NPS-69Gm77081 


tEBRAHT 
TECHNICAL  REPORT  StCTtrm 

'AVAL  POSTGRADUATES,^ 


NAVAL  POSTGRADUATE  SCHOOL 

Monterey,  California 


August  1977 


THE  SECOND-ORDER  SOLUTION  FOR  WAVE 
INTERACTION  WITH  A  SUBMERGED  CYLINDER 
by 
C.  J.   Garrison 
and 
G.  W.  Smith 


Approved  for  public  release;  distribution  unlimited 


FEDDOCS 

D  208.14/2:NPS-69Gm77G81 


NAVAL  POSTGRADUATE  SCHOOL 
Monterey,  California 


Rear  Admiral  I sham  Linder  J.  R.  Borsting 

Superintendent  Provost 


THE  SECOND-ORDER  SOLUTION  FOR  WAVE 
INTERACTION  WITH  A  SUBMERGED  CYLINDER 

The  second-order  solution  of  the  problem  of  the  interaction 
of  a  train  of  regular  waves  with  a  completely  submerged,  hori- 
zontal circular  cylinder  in  water  of  finite  depth  is  presented  for 
two-dimensional  flow.   The  incident  wave  is  developed  as  a  second- 
order  Stokes'  wave  by  use  of  a  perturbation  method  and  the  solution 
of  both  the  first-order  and  second-order  scatter  potentials  is  obtained 
numerically  using  the  Green's  function  approach.   The  hydrodynamic 
pressure  and  resulting  first-order  and  second-order  force  coefficients 
are  determined  numerically  and  presented  for  various  values  of  water 
depth,  cylinder  depth  of  submergence,  and  wave  length. 

The  work  reported  herein  has  been  supported  by  the  National 
Science  Foundation,  Washington,  D.C.   20550. 


SECURITY  CLASSIFICATION  OF  THIS  PAGE  (When  Data  Entered) 


REPORT  DOCUMENTATION  PAGE 


READ  INSTRUCTIONS 
BEFORE  COMPLETING  FORM 


I.     REPORT  NUMBER 

NPS-g9Gm77081 


2.  GOVT   ACCESSION   NO 


3.     RECIPIENT'S  CATALOG  NUMBER 


4.     TITLE  (and  Subtitle) 

THE  SECOND-ORDER  SOLUTION  FOR  WAVE  INTERACTION 
WITH  A  SUBMERGED  CYLINDER 


5.     TYPE  OF   REPORT  &  PERIOD  COVERED 

Final  report;  Sept. -'74 

to  Aug. '77 


6.  PERFORMING  ORG.  REPORT  NUMBER 


7.     AUTHORfaj 

C.J. Garrison 
G.W.Smith 


8.  CONTRACT  OR  GRANT  NUMBERfsJ 

NSF  Grant  No. 
ENG  73-04019  A01 


9-  PERFORMING  ORGANIZATION  NAME  AND  ADDRESS 

Naval  Postgraduate  School 
Monterey,  California  93940 


10.  PROGRAM  ELEMENT,  PROJECT,  TASK 
AREA  &  WORK  UNIT  NUMBERS 


H.  CONTROLLING  OFFICE  NAME  AND  ADDRESS 

Naval  Postgraduate  School 
Monterey,  California  93940 


12.     REPORT   DATE 

1  Aug     1977 


13.     NUMBER  OF  PAGES 
—S3 


14-     MONITORING  AGENCY  NAME  4    ADDRESSf//  dlllerent  from  Controlling  Office) 


15.     SECURITY  CLASS,  (of  this  report) 


unclassified 


I5«.     DECLASSIFICATION/  DOWNGRADING 
SCHEDULE 


16.     DISTRIBUTION   STATEMENT  (of  thla  Report) 


Approved  for  public  release;  distribution  unlimited. 


17.     DISTRIBUTION  STATEMENT  (of  the  abstract  entered  In  Block  30,   If  different  from  Report) 


18.     SUPPLEMENTARY  NOTES 


19.     KEY  WORDS  (Continue  on  reverse  aide  It  necessary  and  Identity  by  block  number) 

wave/body  interaction 
wave  forces 

perturbation  expansion 
Green's  function 


20.     ABSTRACT  I  Continue  on  reverse  aide  It  necessary  and  Identity  by  block  number) 

The  second-order  solution  of  the  problem  of  the  interaction  of  a 
train  of  regular  waves  with  a  completely  submerged,  horizontal  circular 
cylinder  in  water  of  finite  depth  is  presented  for  two-dimensional  flow. 
The  incident  wave  is  developed  as  a  second-order  Stokes'  wave  by  use  of 
a  perturbation  method  and  the  solution  of  both  the  first-order  and 
second-order  scatter  potentials  is  obtained  numerically  using  the  Green's 


DD  ,  ^NRM73  1473 


EDITION  OF   1   NOV  65  IS  OBSOLETE 
S/N    0102-014- 6601   | 


UNCLASSIFIED 


SECURITY  CLASSIFICATION  OF  THIS  PAGE  (When  Data  Sntarad) 


-LCUWlTY   CLASSIFICATION  OF  THIS  PAGECWhwi  D«"  Entered) 


function  approach.   The  hydrodynamic  pressure  and  resulting  first-order 
and  second-order  force  coefficients  are  determined  numerically  and 
presented  for  various  values  of  water  depth,  cylinder  depth  of  submerge: 
and  wave  length. 


UNCLASSIFIED 


SECURITY   CLASSIFICATION  OF  THIS  PAG£fWh»n  Datm  Ei 


ABSTRACT 

The  second-order  solution  of  the  problem  of  the  inter- 
action of  a  train  of  regular  waves  with  a  completely  sub- 
merged, horizontal  circular  cylinder  in  finite  depth  water 
is  presented  for  the  two-dimensional  case.   The  incident 
wave  is  developed  as  a  second-order  Stokes1  wave  by  use  of 
a  perturbation  method  and  the  solution  of  both  the  first- 
order  and  second-order  scattering  potentials  is  obtained 
numerically  using  the  Green's  function  approach.   The  hydro- 
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I.   INTRODUCTION 

The  solution  to  the  two-dimensional  linear  wave/structure 
interaction  problem  has  been  well-studied  for  both  the  finite 
and  infinite  depth  cases.   The  fluid  motion  resulting  from 
the  interaction  of  a  train  of  regular  waves  with  a  submerged 
horizontal  circular  cylinder  in  infinite  depth  water  was 
first  studied  by  Dean  [Ref .  1] .   Ursell  [Ref .  7]  later  studied 
the  problem  anew,  placing  the  solution  on  a  rigorous  basis. 
More  recently,  Ogilvie  [Ref.  6]  reconsidered  the  problem.   He 
computed  the  first-order  oscillatory  forces  and  second-order 
steady-state  forces  for  the  following  cases:   (a)  the  cylinder 
held  fixed,  (b)  the  cylinder  forced  to  oscillate  sinusoidally 
in  still  water,  and  (c)  the  neutrally  buoyant  cylinder  allowed 
to  respond  to  wave  motion. 

It  can  easily  be  shown  that  the  steady-state  part  of  the 
second-order  forces  can  be  computed  from  the  first-order 
potential  alone.   Accordingly,  Ogilvie  did  not  obtain  a 
complete  second-order  solution;  the  steady-state  forces 
arise  from  the  first-order  velocity- squared  terms  in  Ber- 
noulli's equation.   In  addition  to  the  steady-state  forces, 
second-order  oscillatory  forces  also  exist  which  have  a 
frequency  twice  that  of  the  first-order  forces  and  represent 
the  subject  of  the  present  work. 

This  report  reconsiders  again  the  submerged  horizontal 
cylinder  problem,  extending  the  analysis  to  include  water 
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of  finite  depth.   The  cylinder  is  considered  to  be  completely 
submerged  and  held  fixed  in  a  train  of  regular  waves.   The 
objective  is  to  determine  the  complete  second-order  solution 
and,  accordingly,  determine  the  oscillatory  second-order 
forces  as  well  as  the  steady-state  second-order  forces. 

The  problem  is  treated  as  a  regular  perturbation  problem 
in  the  small  parameter,  incident  wave  height/ cylinder  radius. 
Proceeding  in  this  way  the  incident  wave  appears  as  a  second- 
order  Stoke' s  wave.   The  boundary-value  problems  for  both 
the  first-order  and  second-order  potentials  are  established 
and  the  solutions  to  both  are  obtained  by  use  of  the  Green's 
function  method. 
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II.   THEORETICAL  DEVELOPMENT 

A.   FORMULATION  OF  THE  PROBLEM 

A  systematic  representation  of  the  problem  considered  is 
illustrated  in  Fig.  (1) .   A  rigid  right  cylinder  of  radius  a, 
submerged  to  a  depth  d  in  water  of  depth  h,  is  fixed  in 
a  train  of  regular  waves  propagating  in  the  positive  x  direc- 
tion.  The  basic  problem  is  that  of  calculating  the  hydro- 
dynamic  pressure  on  the  immersed  cylinder  and  the  resulting 
force  correct  to  the  second-order  in  the  wave  height  of  the 
incident  wave. 

Assuming  the  fluid  to  be  irrotational ,  a  velocity 
potential,  §,    may  be  defined  as: 

q  =  V$(x,y,t)  (1) 

where  q  denotes  the  fluid  velocity  vector.   (The  barred 
quantities  denote  dimensional  quantities.)   Moreover,  assuming 
the  fluid  to  be  incompressible,  and  homogeneous,  the  velocity 
potential  must  satisfy  the  Laplace  equation, 


V2$(x,y,t)  =  0  (2) 


within  the  fluid  region. 

In  addition  to  the  differential  equation,  Eq.  (2) ,  $ 
must  satisfy  certain  boundary  conditions.   In  specific, 
these  are: 
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1.  The  velocity  of  the  fluid  normal  to  the  bottom 
boundary  must  vanish. 

2.  The  velocity  normal  to  the  surface  of  the  cylinder, 
defined  by  S,  (x,y)=0,  must  vanish. 

3.  The  kinematic  boundary  condition  on  the  free  surface, 
defined  by  S~ (x,y , t) =n (x, t) =0 ,  must  be  satisfied. 

4.  The  dynamic  boundary  condition  on  the  free  surface 
must  be  satisfied. 

These  boundary  conditions  may  be  stated  mathematically  as  follows 


>-(x,-h,t)=0  (3) 


V$«VS,=0  (4) 


3S 
q-S2+   ^=0  (5) 


>-(x,n,t)+  *5[>-(x,n,t)  +  $-(x,n,t)  ]+gn(x,t)=gK  (6) 


where  n  denotes  the  elevation  of  the  free  surface,  g  denotes  the 
acceleration  of  gravity,  and  K  denotes  the  Bernoulli  constant. 
For  convenience  in  carrying  out  the  solution,  the  variables 
are  next  made  dimensionless  using  the  cylinder  radius,  and  the 
wave  frequency,  a  ,  as  follows: 
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x  =  x/a      y  =  y/a      d  =  d/a      h  =  h/a 


H  =  H/2a     K  =  K/a      v  =  a2a/g    t  =  at 


<J>  =  a$/ga    ri  =  n/a 


(7) 


H  denotes  the  dimensionless  wave  height,  K  denotes  the 
dimensionless  Bernoulli  constant,  and  t  denotes  the 
dimensionless  time. 

Utilizing  the  parameters  defined  in  Eq.  (7),  the  boundary 
value  problem  defined  in  Eqs.  (2-6)  may  be  rewritten  concisely 
in  dimensionless  form  as: 


V  <J>(x,y,t)  =0       in  the  fluid  region 


(8) 


<J>  (x,-h,t)  =  0 


(9) 


<J>n(x,y,t)  =0        on  S1(x,y)  =  0 


(10) 


<frx(x,n,t)nx(x,t)  -  <j>  (x,n,t)  +  vnt(x,t)  =  o 


(ID 


<t>t(x,n,t)  +  jUx(x,n,t)2  +  <$>y(x,r},t)2]   +  n(x,t)  =  k  (12) 


B.   PERTURBATION  EXPANSION 

Expanding  <j>  and  n  in  terms  of  a  perturbation  parameter, 
z ,  provides  a  means  for  converting  the  nonlinear  boundary-value 
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problem  into  a  series  of  linear  problems.  The  solution  to 
each  linear  problem  may  be  tractable  whereas  the  nonlinear 
problem  is  unsolvable. 

Since  <J>,  r\ ,  and  K  are  functions  of  the  small  parameter, 
e,  they  may  be  written  in  the  power  series: 


<Mx,y,t;e)  =  £  en6  (x,y, t)  (13) 

n=l    n 


n(x,t;e)  =  £  enn  (x,t)  (14) 

n=l    n 


and 


K=£?:nKn  (15) 

n=2 


In  Eqs.  (11)  and  (12),  <j>  contains  r],    and,  therefore,  e., 
implicitly.   This  can  be  converted  to  an  explicit  form  in 
e  by  use  of  the  Taylor  series  expansion: 


0(x,n,t)  =  £  n(x't;£)m  C3m^(x,y,t)3v-n         (16) 
m=0    ml      3y         y  J 


The  perturbation  parameter,  e,  is  related  to  the  wave  height 
in  an,  as  yet,  unknown  manner. 

Equations  (13-16)  as  well  as  appropriate  derivatives 
similar  to  Eqs.  (13-16)  may  now  be  substituted  into  the 
boundary-value  problem  given  in  Eqs.  (8-12)  to  obtain  a 
series  of  linear  boundary-value  problems  for  <f>.  ,  <f>_,  etc. 


20 


That  is,  upon  substitution  and  equating  coefficients  of  like 
powers  of  e ,  the  following  problems  for  the  first  two  terms 
in  the  expansion  for  4>  can  be  precipitated: 
First-order  boundary-value  problem  (e): 


V2cj)1(x/y/t)    =    0  (17) 


<\>1    (x,-h,t)    =    0  (18) 

*ln(x'y,t)    =    °  °n   Sl(x'y)    =    °  (19) 


*1  (x,o,t)   -  vnlt(x,t)   =  o 


V24)2(x,y,t)    =    0 


(20) 


nx(x,t)    +    4»lt(x,0,t)    =   0  (21) 


2 

Second-order   boundary-value   problem    (e    ): 


(22) 


<J)2y(x,-h,t)    =    0  (23) 

<f>2n(x,y,t)    =0  on   S1(x,y)    =   0  (24) 

<t>2    (x,o,t)    -  vn2t(x,t)   = 

(25) 

-n1(x,t)<|>1     (x,o,t)   +  <J>lx(x,o/t)nlx(x,t) 
n2(x,t)   +  $2t(x,o,t)   =  -n1(x,t)<j>1  t(x,o,t)    - 

nlt(x,t)<|>ly(x,0,t)     -   ~[<Plx(xf0ft)2    +    <J>ly(x,0,t)2]    +   K2     (26) 
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There  is  a  striking  similarity  between  the  first-order 
and  second-order  problems;  only  the  right  hand  side  of  the 
free  surface  boundary  conditions  differ.   In  the  first-order 
problem  the  free  surface  boundary  conditions  are  homogeneous 
whereas  the  second-order  free  surface  boundary  conditions 
are  non-homogeneous,  the  right  hand  side  being  functions  of 
the  first-order  velocity  potential,  first-order  free  surface 
elevation,  and  their  derivatives. 

The  first-order  potential  may  be  represented  by  a  function 
which  is  periodic  with  the  fundamental  frequency.  Accordingly, 
the  complex  potential,  u,  (x,y)  ,  may  be  expressed  as: 


♦1(x,y#t)  =  abRe  [iu1(x,y)e  lfc]  (27) 


where  R   denote  the  real  part,  a  =  2ua/L,L  denoting  the  wave 
length,  and  b  is  an  unknown  real  constant. 


C.   THE  FIRST-ORDER  BOUNDARY  VALUE  PROBLEM 

Since  the  first-order  boundary-value  problem  is  linear,  the 
total  first-order  potential  may  be  expressed  as  the  sum, 

}1   -  <f>*  +  <j,S  (28) 

I  .  S 

where  <f>,   denotes  the  incident  wave  potential  and  <J>,   denotes 

the  scatter  potential  which  is  due  to  the  presence  of  the  cylinde 

The  space  dependent  part  of  the  complex  potentials  are  expressed, 

in  view  of  Eqs .  (27)  and  (28)  as, 
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u1  =  u*  +  u^  (29) 


The  potential,  u,  ,  represents  only  a  regular  wave  propa- 
gating along  the  x-axis  and  this  potential  must  satisfy  the 
first-order  boundary-value  problem  when  no  rigid  body  is 
present.   Therefore,  <J>,  satisfying  Eqs .  (17-18)  and  (20-21), 
is  given  by 


I    _  1  cosh  [a(y+h)]  eiax 
1      a     cosh  (ah) 


The  parameters  a  and  v  are  related  through  the  well-known  equation 
from  linear  wave  theory. 

2- 

v  =  2-£-  =  a  tanh(ah)  (31) 

g 

which,  in  dimensional  form,  using  a  =  ka ,  and  k  =  2tt/L  is: 

a2=  gk  tank(kh)  (32) 

Since  the  solution  for  u.  is  known,  the  boundary-value 
problem  for  u,  may  now  be  established.   Substituting  Eqs.  (30) , 
(29)  and  (27)  into  the  first-order  problem  given  by  Eqs.  (17-21) , 
and  eliminating   n-,  between  Eqs.  (20)  and  (21)  results  in  a 
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boundary-value  problem  for  the  first-order  scatter  potential 
as  follows: 


V2u^(x,y)  =  0 


(33) 


u*  (x,-h)  =  0 


(34) 


S  ,    .        1 

uln(x,y)  = 


cosh  (ah) 


n   sinh[a(y+h)]  +  in   cosh[a(y+h) 

.  y  x 


iax 


on  S.(x,y)  =  0 


(35) 


u^y(x,0)  -  vu^(x,0)  =  0 


(36) 


where  n   and  n   denote  the  x-  and  y-components ,  respectively, 
of  the  unit  normal  vector,  n  =  in   +  jn  ,  directed  outward 
into  the  fluid. 

D.   SECOND-ORDER  BOUNDARY-VALUE  PROBLEM 

In  view  of  the  linearity  of  the  second-order  problem, 
the  same  procedure  as  used  in  the  first-order  problem  may 
be  applied  and  leads  to  the  representation  of  the  second-order 
potential  as  the  sum: 


<f>2    -    <j>2   +    <j>2 


(37) 


where  <j>2  denotes  the  second-order  incident  wave  potential 

S 
and  (J> 2  denotes  the  scatter  potential. 


24 


Considering  the  case  where  there  is  no  body  present, 

S     S 
there  is  no  scattered  wave  and  therefore,  <f>,  =  ^  =  0 . 

Additionally,  the  boundary  conditions  on  the  immersed  cylinder, 
Eqs.  (19)  and  (24),  are  not  applicable.  Thus,  when  the  substitu- 
tions of  4>.  for  <£>,  and  <j>2  for  cj>2  are  made  into  Eqs.  (22)  , 
(23),  (25)  and  (26),  along  with  Eqs.  (20),  (21),  (27)  and 
(31),  and  upon  eliminating  n,  and  ru  between  Eqs.  (25)  and 
(26) ,  the  resulting  boundary-value  problem  for  the  second- 
order  incident  wave  potential  becomes: 


V24>2(x,y,t)  =  0  (38) 


<J>2  (x,-h,t)  =  0  (39) 

4>2  (x,0,t)  +  v<J>2tt(x,0,t)  =  -  |b2(a2-v2)Re[iei2(ax~t)  ] 

(40) 

The  periodic  solution  to  the  incident-wave  boundary-value 
problem  specified  by  Eqs.  (38-40)  is  known  and  may  be  expressed 
as 


<J>2    =  |ab2vRe[iu2(x,y)e"l2t]  (41) 


where  the  second-order  complex  potential,  u? ,  is  given  by: 


I(x,y)    =   .     1    cosh[2a(y+h)l    ei2ax  (42) 

sinh    (ah) 
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Eqs.  (41-42)  are  familiar  in  wave  theory  and  are  exactly  the 
same  form  as  that  given  by  Ref.  4  for  second-order  Stokes' 
waves. 

Having  developed  a  solution  for  the  second-order  incident 

wave  potential,  $2,  the  problem  for  second-order  scatter 

I     S 
potential  may  be  determined.  To  this  end,<J>_=$,  +  <j>..  and 

I  S 
*2=<^?  +  ^2  are  sut>stituted  into  Eqs.  (22-26),  and  the  pre- 
viously determined  solutions  for  the  incident  wave  potentials, 
Eqs.  (30)  and  (42),  are  utilized.  After  applying  Eqs.  (20) , (21) , 
(27),  (29)  and  (41)  and  eliminating  r\  between  Eqs.  (25)  and  (26) 
the  resulting  boundary-value  problem  for  the  second-order  scatter 
potential  becomes: 


V2<f>2(x,y,t)    =    0 


(43) 


<J>2    (x,-h,t)    =    0 


(44) 


4>2n(x,y,t) 


ab2    R 


sinh    (ah) 


|(n      cosh[2a(y+h) ] 


-    in      sinh [2a (y+h) ] ) e 


i2 (ax-t) 


(45) 


on  S. (x,y)    =   0 


*2y  +    *2tt 


a2b2    D    r .  ,1  S      S  IS  .   I      S  ,.-. 

— = —  R    [i(— u,    u,         +   u,u,         -    6u,    u,  (46) 

2         e        v   ly   lyy  1   lyy  ly   ly         K      ' 


1SI  .    I   ,  S  0S2         -.s2a    -i2t1 

+  — u,    u,         -   4u,    u,      -   2u.         -    3u,       ) e  ] 

v   ly   lyy  lx   lx  lx  ly  J 


+  U(x) 


on   y   =   0 
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where  U(x)  is  a  non-periodic  function  generated  by  the 
substitution  of  <j>,  into  Egs .  (2.5)  and  (26).   For  brevity, 
this  function  is  not  written  out  since  it  is  not  needed  in 

determining  the  second-order  pressures  and  forces. 

S 
The  boundary-value  problem  in  <j>~  given  in  Eqs,  (43-46) 

is  time  dependent  according  to  e     ,  except  for  the  U(x) 

c 

term  in  Eq.  (46)  .   Therefore,  the  solution  for  <|>.  may  be 
taken  in  the  form 

♦  5  -  |  ab2vRe[i[u^  (x,y)e~l2t  +  u|(x,y)]J     (47) 

where  the  last  term  denotes  the  time-independent  portion  of 

the  complex  potential.   Separate  boundary-value  problems  for 

S       "  S 
u?   and  u?  arise  from  the  substitution  of  Eq.  (47)  into  Eqs. 

S 
(43-46) .   However,  as  will  be  demonstrated  only  the  $~.     term 

is  required  for  determination  of  the  pressure  to  the  second 
order.   Therefore,  the  time-derivative  of  the  time-independent 
part  of  <j>  ?   will  be  zero,  and,  accordingly,  will  not  be  developed 
further  here;  the  solution  for  u2  only  will  be  pursued. 

When  Eq.  (47)  is  substituted  into  Eqs.  (43-46),  the  result- 
ing boundary-value  problem  for  the  second-order  scatter  potential 
takes  the  form: 

V2u^(x,y)  =  0  (48) 


u^y  (x,-h)=  0  (49) 
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s 
u-    (x,y)    = 


/n  sinli     (ah)     L 


n      sinh[2a(y+h) ] 

y 


+   inx  cosh[2a(y+h)  ]1    el2ax     on   S1(x/y)=0        (50) 


u2y(x,0)    -    4vu2(x,0)    =    f*(x)  (51) 


where 


.„,,    .         2arl    S      S  IS  ,   I      S  IS      I  ,,- 

f*  (x)    =    3V^ulyulyy  +   u^^y   -    6ulyuly  +  ^ulyulyy  (52) 

2  2 

,    I      S  -    S  ,    S      , 

-    4u,     u,        -    2lln  -    3u.        J„_n 

lx  lx  lx  ly     y-0 


There  is  considerable  similarity  in  the  first-order 
and  second-order  problems,  the  only  difference  in  form  being 
that  the  second-order  free  surface  boundary  condition, 
Eq.  (51),  is  non-homogeneous. 

Because  of  the  non-homogenity  of  Eq.  (51),  further  use  is 
made  of  linear  superposition  defining  u»  as  the  sum, 


u2  =  u2   +  u2  (5  3) 

S°  S* 

where  u,,   denotes  the  homogeneous  solution  and  u2   the 

particular  solution  to  the  boundary-value  problem  as  stated 

S*      S° 
in  Eqs.  (48-51).   More  precisely,  u2   and  u2   are  defined  as 

the  solutions  to  the  boundary-value  problems  obtained  from 

the  substitution  of  Eq.  (53)  into  Eqs.  (48-51)  as  follows: 
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s  * 
Particular  Solution  (u^ 


V2U2*(x,y)  =  0  (54) 


u^*  (x,-h)  =  0  (55) 


U2*(x,0)  -  4vu2*(x,0)  =  f*(x)  (56) 


S° 
Homogeneous  Solution  (u_  ) 


V2uf(*,y)    =   0  <57> 

s° 

Uj    (x,-h)    =   0  (58) 


S°  S° 

u^    (x,0)    -    4vu^    (x,0)    =    0  (59) 


s°  i  r 

u      (x,y)    =   j nv   sinh[2a(y+h)  ] 

sinh  (ah)   L  y 


* 


+  inx  cosh[2a(y+h)  ]lel2ax  -  ^(x^) 


(60) 


on  S- (x,y)  =  0 

By  the  division  of  the  second-order  problem  into  homogeneous 

S* 
and  particular  solutions,  the  non-homogeneous  problem  for  u« 

contains  no  boundary  condition  on  the  cylinder  surface.   The  result- 

S* 
ing  boundary-value  problem  for  u2  ,  as  stated  in  Eqs.  (54-56),  is 

identical  to  that  associated  with  the  linear  problem  for  the  potentia 
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resulting  from  a  harmonic  pressure  variation  of  amplitude 

distribution  f*(x)  on  the  free  surface  in  water  of  depth  h. 

S° 

The  homogeneous  boundary-value  problem  for  u„  ,  defined  by 

Eqs.  (57-60)  ,  is  similar  in  form  to  the  first-order  problem 
given  by  Eqs.  (33-36),  the  significant  difference  being  the 
term  4v  in  place  of  v  which  occurs  in  the  free  surface  boun- 
dary  condition.   Thus,  the  method  of  solution  for  u2   will 
be  similar  to  that  of  the  first-order  scatter  potential 


E.   DEFINITION  OF  THE  INCIDENT  WAVE 

Having  developed  the  first-order  and  second-order  boun- 
dary-value problems,  it  is  now  appropriate  to  completely 
define  the  incident  wave  height  and  in  the  process  determine 
expressions  for  the  unknown  constants  b  and  K_ .   Solving 
Eqs.  (21)  and  (26)  for  n,  and  n2'  respectively,  and  then 
substituting  the  results  into  the  expression  for  the  free 
surface  elevation  as  given  by  Eq.  (14)  yields: 


n(x,t)  =  -  E*lt  +  €2[-*2t  +  *lt*lyt  +  *ltt*ly       (6D 


I(*lx  +  *ly>  +  K2!  +  0(e3) 


where  the  potentials  are  evaluated  at  y  =  0 .   Eq.  (61) 
expresses  the  free  surface  elevation  in  terms  of  the  total 
potential,  and  therefore,  includes  the  effects  of  both  the 
incident  and  the  scattered  waves.    Hence,  as  it  is  of  inter- 
est to  obtain  an  expression  for  the  free  surface  elevation 
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of  the  incident  wave,  the  scatter  potentials  are  set  to 

S     S 
zero,  i.e.  <f>_  =  <j>   =  0.   Evaluating  Eq.  (61)  using  the  known 

solutions  for  the  incident  wave  potentials,  Eqs.  (30)  and 

(42),  along  with  Eqs.  (27)  and  (41),  then  yields: 


H^t)  =  ebRje1^"^]  +  s2[b2(v^a2)  +  K2         (62) 

ab  cosh  (ah)  [2  +  cosh  (2ah)  ]      i2 (ax-t) , 
+  -2 = R  le        ] 

sinhi  (ah)  e 

+  0(e3) 


where  n  (x,t)  denotes  the  dimensionless  surface  elevation 
for  the  incident  wave  with  no  cylinder  present.   If  the  x-axis 
is  placed  at  the  mean  water  line  then  the  constant  term  must 
vanish  and,  therefore, 


u2/  2     2. 
m   b  (a   -  v  )  (63) 

2  4  v 


The  remaining  terms  in  Eq.  (62)  are  time -dependent,  periodic 
functions,  the  second-order  term  having  a  frequency  of  twice 
that  of  the  first-order. 

Defining  the  dimensionless  wave  height  given  in  Eq.  (7) 
as  the  difference  in  elevation  between  the  crest  and  trough 
of  the  incident  wave,  eb  must  represent  the  wave  amplitude. 
The  second-order  term  in  Eq.  (62)  having  frequency  twice  the 
fundamental  frequency  makes  equal  contributions  to  surface 
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elevation  at  both  the  crest  and  trough,  so  that  it  contributes 
nothing  to  the  wave  height.   Thus,  in  terms  of  dimensionless 
wave  height,  b  is  given  by 

H  =  eb  =  H/2a  (64) 

where  H  denotes  the  elevation  difference  between  the  wave 
crest  and  trough. 

Expressing  the  incident  wave  profile  in  terms  of  the 
dimensionless  wave  height  yields: 


nI(x,t)  =  H  cos  (ax-t)  (65) 


2 

H  a  cosh(ah)  r~  ,     ,  /,v,»  ,      r„  ,    ,  .  , 
+ [2  +  cosh  (ah)  ]  cos  [2  (ax-t)] 

4  sinh  (ah) 


It  may  be  noted  that  Eq.  (65)  agrees  with  the  expression 
for  the  second-order  Stokes1  wave  given,  for  example,  by 
Ref.  4. 

F.   PRESSURES  AND  FORCES 

The  pressure  may  be  formulated  by  use  of  Bernoulli's 
equation  arranged  as  follows: 


2      2 

P(x,y,t)  =  -p$r-  -  ip[«-  +  «-  ]  -  pgy  +  pgK  (66) 

*"       ^       X  V 


where  P(x,y,t)  denotes  the  dimensional  pressure  and  p  denotes 
the  fluid  density.   By  use  of  Eqs.  (7),  (13-16),  (27),  (41), 
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(47) ,  and  (64) ,  Eq.  (66)  may  be  reduced  to  dimensionless  form, 
and  carried  to  the  second-order  in  wave  height.   The  result  is 


p(x,y,t)  =  -y  -  HaR^u^"1*1]  (67) 


2  2 


Ha 


r6\T        2      2,  -i2t, 
R^  t-— — un  -  u,   -  u,  )e     ] 


4v  L  e   a   2    "lx   "ly' 


2        2     2-. 

+  Iulxl   +  |uly|   +-2--  lj 


In  Eq.  (67)  p(x,y,t)  denotes  the  dimensionless  pressure 
coefficient  and  is  defined  as: 


p(x,y,t>  =  L<*-£i*J  (68) 

pga 


The  first  term  in  Eq.  (67)  represents  the  hydrostatic 
pressure  as  y  is  the  dimensionless  depth  beneath  the  mean 
free  surface  (y  =  0) .   The  second  and  third  terms  are  harmonic, 
the  third  term  having  twice  the  frequency  of  the  second,  and 
representing  the  harmonic  second-order  contribution.   The 
remaining  terms  in  Eq.  (67)  are  independent  of  time  and  result 
in  the  time-average  or  steady-state  force  components. 

Expressing  the  components  of  the  wave  force  vectors  in 
terms  of  integrals  of  the  pressure  over  the  cylinder  surface 
area,  the  dimensionless  force  coefficients  may  be  written 
as 


^(t)  =  -  yp(x,y,t)ni  dS1,   i  =  1,2  (69) 

Sl 
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where  the  dimensionless  force  coefficients  in  the  x  and  y 
directions  are  defined,  respectively,  as 


F  (t) 
Cx(t)  =  -2-g-  (70) 

pga 


F  (t) 
C2(t)  =  -£jg-  (71) 

Pga 


where  F  (t)  and  F  (t)  denote  the  force  components.  Addi- 
x         y 

tianally,  dS-,  denotes  a  dimensionless  differential  arc 
length  along  the  circumference  of  the  cylinder. 
Applying  Eq.  (67)  to  Eq.  (69)  yields: 


Ci(t)  =  -  y   /n±dS1  -  Ha  f^i^e   lt]n±  dS1 


(72) 


Sl  Sl 


„2ra2         /I     r6v2  2  2.     -i2t 

"   H    [4v"        /Re[-T"u2    -   ulx  "   uly)e 
bl 


2         2     2 

+  |ulx'   +  |uly'   +  ^T  '  l]nidS1  +  0(H3) 

a 


The  first  term  in  Eq{72)  results  from  the  hydrostatic 
pressure  on  the  cylinder  so  it  represents  simply  the  buoyancy 
force.   Disregarding  the  hydrostatic  pressure  the  hydrodynamic 
force  may  be  written  as 
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Ci(t)  =   H  F1±    cos(61:L-t)  +  H2  (73) 


+  H2[F2i  cos(S2i-2t)  +  F^?]  +  0 (H3) 


where  the  first-order,  second-order,  and  steady-state  force 
coefficients  and  phase  shift  angles  are  defined  by  comparison 
of  Eqs.  (72)  and  (73)  as  follows: 


\ie  =    a  e  7UinidSl 


i<5. 

F, .e   "  =  a   /u.n.dS,  (74) 

S. 


1 

S 


F2iel  ^  =  h       /(£i"u2  "  V'  -  O^dS,  (75) 


2      2 

lx  "  ulyyiii~~l 
1 


F2i  =  ^q/(lulx|2  +  l^lyl2  +  4-  »^1  <76> 


The  dimensionless  force  coefficients,  F. .  and  F-.  are  real. 

The  first  term  in  Eq.  (73)  represents  the  first-order 
forces  of  the  fundamental  frequency,  a.   The  periodic  portion 
of  the  second  term  in  Eq.  (73)  represents  the  second-order 
contribution  to  the  force  at  twice  the  fundamental  frequency. 
The  last  second-order  term  represents  the  steady-state  or 
time- independent  contribution,  sometimes  called  drift-force. 

Evaluation  of  the  force  coefficients  defined  by  Eqs. 
(74-76)  is  the  main  objective  of  this  work.   Therefore, 
it  is  necessary  to  evaluate  the  first-order  and  second-order 
complex  potentials,  u,  and  u2 ,  as  well  as  derivatives  of  u. 
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on  the  surface  of  the  cylinder.   Additionally,  evaluation 
of  u,  and  its  derivatives  on  the  mean  free  surface  (y  =  0) 
is  required  in  order  to  carry  out  the  solutions. 
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III.   METHOD  OF  SOLUTION  AND  NUMERICAL  PROCEDURES 

A.   SOLUTION  OF  THE  FIRST-ORDER  PROBLEM 

The  use  of  a  Green's  function  to  express  the  solution 
to  the  first-order  boundary-value  problem  was  first  formu- 
lated by  John  [Ref .  5] ,  and  applied  to  submerged  ellipsoids 
by  Garrison  and  Rao  [Ref.  3].   This  method  is  considered  to 
be  applicable,  in  principle,  to  arbitrary  shapes  and  is 
mathematically  the  most  straightforward.   The  Green's  func- 
tion, G,  of  unit  strength  which  satisfies  Eq.  (33)  as  well 
as  the  boundary  conditions,  Eqs .  (34)  and  (36)  is  given  by: 


G(x,y;£/n;v)  =  In  r  -  In  r'  +  2PV  f 


00  _ 


cosh [y (y+h) ] cosh [y (n+h) 3 
_(coshyh) (vcoshyh-ysinh  yhT 


-  e 


-yh  sinh(un)  smh(uy)  1     i   rij 

■  r1 1    u  \ —      cos  x-E  dy 

smn  yh      J     i   •■  i  p 


(77) 


4tt 


1  2a1h+sinh(2a1h)  COSh [al (y+h) ] COSh [al (n+h) ] COS [al ' X"C I ] 


where : 


r  -  [(x-o2   +  (y-n)2]2 


(78) 


r'  =  [(x-c)2  +  (y+n)2]2 


(79) 


The  symbol,  a..,  is  defined  in  terms  of  h  and  v  as  the 
solution  to  the  equation 


37 


F(a1,h,v)  =  0  (80) 


where  F  is  defined  by 


F(a.,h,v)  =  a  tanh(a.h)  -  v  (81) 


Comparison  of  Eqs .  (31),  (80)  and  (81)  demonstrates  that  a, 
is  clearly  the  equivalent  of  a.   However,  this  will  not  be 
the  case  for  the  second-order  problem,  thus,  requiring  the  use 
of  separate  notation.   In  Eq .  (77),  PV  denotes  the  principal 
value  of  the  integral. 

The  Green's  function  given  in  Eq .  (77)  was  also  given  in 
series  form  by  John  [Ref.5]  as 

2   2 
«  (vk+v  )cos[yk(y+h) ]  -y .  | x-£ 

G(x,y;^,ri;v)  =  2tt  2^  75 2 cos  [yk (n+h) ] e 

k=l  vPk-hVk(lk+v  ) 

iajx-sl 
47TCOsh  [a.  (y+h)  ]  cosh  [a,  (n+h)  ]e 

-  i = (82) 

2a  h  +  sinh  (2axh) 

where  a^  is  as  defined  in  Eq.  (80)  and  the  quantities  y. 
are  defined  as  the  real  positive  roots  of  the  equation 


K(yk,h,v)  =  0  (83) 


where  the  function  K  is  defined  as 
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K(y,h,v)  =  y  tan(yh)  +  v  (84) 


S 
Following  the  Green's  function  method  of  solution,  u, 

is  written  as  the  integral  over  the  cylinder  arc  length, 

S  ,  as: 


wjfi^>* 


U1(X/Y)  =  ItF   /fiUrn)G(xfy;5,n;v)  dSx      (85) 
Sl 

where  (£,n)  denotes  points  on  the  immersed  surface,  f.(£,n) 
denotes  the  unknown  source  strength,  and  dS ,  =  ds",  /a  denotes 
the  differential  arc  length  on  the  cylinder  surface,  made 
dimensionless  with  the  characteristic  length,  a. 

Although  the  solution  to  the  first-order  boundary-value 
problem  as  stated  in  Eqs.  (33-36)  is  given  by  Eq.  (85) ,  the 
source  strength  function,  f.(£,n),  must  be  determined  in 

order  to  evaluate  the  potential.   From  potential  theory, 

S 
the  normal  derivative  of  the  potential,  u, (x,y) ,  on  the 

surface  of  the  cylinder  is 

uln(x'Y)  =  1   fi(X'y)  +  J?   /f1U,n)Gn(x,y;£,n7v)dS1  (86) 

Sl 

where  G   (the  normal  derivative  of  G)  may  be  determined  by 
differentiation  of  either  Eq.  (77)  or  (82)  in  a  straight- 
forward manner . 

Applying  the  final  boundary  condition,  Eq .  (35),  leads 
to  an  integral  equation  which  may  be  solved  for  f. , 
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j  f1(x,y)   +  ~      /f1(C,n)Gn(x,y?^n,v)    dsx 


(87) 
j 

si 


cos 


_^y|n   sinh[a(y+h)  ]    +   in^osh  [a  (y+h)  ] 


lax 
e 


The  solution  for  f,  from  Eq .  (87)  may  then  be  used  in  Eq. 
(85)  to  determine  the  potential,  u, . 

B.   SOLUTION  OF  THE  SECOND-ORDER  PROBLEM 

The  similarity  in  form  of  the  boundary-value  problem 
for  the  homogeneous  part  of  the  second-order  potential, 
Eqs.  (57-60) ,  to  the  first-order  potential  is  now  utilized 

Since  the  only  significant  difference  occurs  in  Eq.  (59) , 

s° 

where  v  is  replaced  by  4v,  we  may  represent  u9   in  a  form 

similar  to  Eq.  (85)  as 


s°.   _.  _  i     L 

2tt  J 


U2  (x,Y)  =  fi        /f2U,n)G(x,y;£,v;4v)  dS1  (88) 


Sl 


where  G  (x,y;  E,,r\ ;  4v)  is  defined  by  replacing  v  by  4v  in  Eqs. 
(77)  and  (82).   Therefore,  a2  is  defined  by 


F(a2,h,4v)  =  0  (89) 


and  a^,  is  replaced  by  a2  in  Eqs.  (77)  and  (82), also.   In  the 
case  of  Eq.  (82) ,  y.  is  defined  as  the  real  positive  roots  of 
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K(yk/h,4v)  =  0  (90) 


The  source  strength  function,  f2(£,n)/  appearing  in 
Eq.  (88)  is  again  determined  by  applying  the  kinematic 
boundary  condition  on  the  cylinder  surface  as  given  by  Eq . 
(59) .   The  integral  equation  for  f2  may  then  be  given  by 

1   f2(x'y)  +  T?   yf2(C'n)Gn(x'y;C'n;4v)  dSl  (91) 


Sl 


n    sinh[2a(y+h) ]    +    in   cosh [2a (y+h) ] 

y   5c  


i2ax    S*     , 
e     -  u2n(x'y) 


smh  (ah) 


on  S1(x/y)  =  0 


S* 
However,  to  solve  Eq.  (91)  for  f 2 ,  first  a  solution  for  u2 

S* 

must  be  obtained  and  u~   evaluated  on  the  cylinder  surface. 

2n 

S* 
Thus,  at  this  point  the  solution  to  u2   is  sought. 

S* 

The  complex  potential  u2   is  defined  as  the  solution  to 

the  boundary-value  problem,  Eqs .  (54-56) .   The  form  of  this 
problem  is  recognized  as  being  equivalent  to  that  of  fluid 
motion  produced  by  a  periodic  pressure  distribution  (as 
specified  by  f*(£))  with  frequency  2a   on  the  free  surface.   No 
cylinder  is  considered  to  exist  in  the  fluid  region. 

The  solution  to  this  problem  may  be  formulated  as  an 
integral  over  the  free  surface  extending  from  negative 
infinity  to  positive  infinity.   Using  the  present 
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s* 

dimensionless  representation,  u2   becomes: 

u2*(x,y)  -  ~p   /f*(C)G*(x,y;C/0;4v)  d£  (92) 

2 

The  Green's  function,  G* , for  this  problem  is  given  in  Ref.  8. 
Upon  transforming  the  solution  to  the  non-dimensional  form 


G*(x,y;5,0;4v)  =  -2PV   /"  c°sh[u (y+h) ]cos [y (x-g) ]dp 

-co  4vcosh(i_ih)    -    ysinh(yh) 


47Tcosh  [a9h]  cosh  [a0  (y+h)  ]  cos  [a9  (x-£)  ] 

-i  1 1 i (93) 

2a2h  +  sinh(2a2h) 


where  a2  is  defined  by 


F(a2,h,4v)  =  0  (94) 


S* 

The  normal  derivative  of  u9   is  required  in  order  to 

completely  define  the  kinematic  boundary  condition  on  the 
cylinder  as  given  in  Eq.  (60) .   From  Eq .  (92),  this  becomes 

u2*(x,y)  =  ~   /f*(OG*(x,y;£,0;4v)  d£       (95) 

H 

The  derivative  of  G*  in  the  direction  normal  to  the  cylinder 
surface  may  be  obtained  by  differentiation  of  Eq .  (9  3)  in  a 
straightforward  manner.   Thus^ using  the  values  for  f*(£), 
as  defined  in  terms  of  the  first-order  solution  only, and  as 
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s*      s* 
given  in  Eqs .  (52)  and  (56),  u2   and  u2   may  be  determined 

S*     S° 

Using  u2n  /  u2   may  then  be  determined  to  complete  the 

second-order  boundary-value  problem  solution. 


C.   NUMERICAL  METHODS 

Determination  of  the  forces  on  the  cylinder  surface 
requires  solving  for  the  potentials  u,  ,  u2  and  the  derivatives 
of  u,  at  points  on  the  surface,  as  shown  in  Eq.  (72).   Thus, 
the  problem  is  now  one  of  solving  for  both  the  first-order 
and  the  second-order  scatter  potentials,  and  requires  solving 
for  the  source  strength  functions  f,  and  f„,  as  well  as  the 
free  surface  pressure  distribution  source  strength  function, 
f  .   In  view  of  the  complexity  of  the  equations  it  is  natural 

to  attempt  a  numerical  solution. 

S 
The  first-order  scattering  potential,  u, ,  as  given  by 

Eq.  (85)  is  dependent  upon  the  first-order  source  strength 

function,  f, .   Thus,  it  is  necessary  to  solve  Eq.  (87)  for 

f,  to  determine  u..  .   A  numerical  solution  may  be  developed 

by  dividing  the  cylinder  surface  into  elements  of  length 

A9  =  2ir/m,  with  the  center  of  each  element  assigned  an  index. 

Since  f, (£,n)  is  recognized  as  a  well-behaved  function  for  a 

smooth  surface  cylinder,  it  is  appropriate  to  define  parameters 


a..(v)  =-     /G  (x. ,y. ;£,n;v)  dS  (96) 

ID       tt  Ae.  J    n   i   i  1 


ASlj 


i, j  =  1 ,2  ,3,  ...  ,m 
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hi  ■  cosh(ah)   [ny(xi'yi'  "^ta(yi+h)]  (97) 


+  in  (x.,y.)  cosh[a (y .+h) ] 

X    JL     JL  -J- 


iax. 
e    x 


fu-  ^(x.^)  (98) 

and, thus, Eq.  (87)  may  be  approximated  by  the  complex  matrix 
equation, 


fli  +  aij(v)flj  =  2hlj     i/j  =  1'2""'m  (99) 


Upon  evaluation  of  a. . (v)  using  Eq .  (96),  the  inversion  of 
Eq.  (99)  may  be  carried  out  on  the  digital  computer  to 
determine  f, .  at  each  nodal  point  on  the  surface  of  the 
cylinder. 

For  the  purpose  of  evaluating  a  and  $,  either  Eq .  (77) 
or  Eq.  (82)  may  be  used.   For  evaluation  of  a  and  6  when  i  =  j 
Eq.  (77)  must  be  used  to  allow  separate  treatment  of  the 
logarithmic  singularity  as  r  approaches  0.   The  difficulty 
encountered  with  the  logarithmic  singularity  may  be  overcome 
by  carrying  out  its  integration  analytically.   Garrison 
[Ref.  2]  showed  for  the  calculation  of  a. .  that  the  contribu- 
tion of  the  normal  derivative  of  the  In  (r)  in  Eq.  (77) 
integrated  over  the  singular  element  of  A0  is  A0/2,  not  only 
for  the  nodal  point  (x.  ,y.)  =  (£,n) /  but  for  all  nodal  points 
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(x.,y.)  on  the  cylinder  surface.  Additionally,  to  calculate 
8. .  Garrison  has  developed  a  numerical  approximation  for  the 
integral  of  the  In  (r)  over  the  singular  element  of  A6. 

A  second  singularity  arises  in  the  evaluation  of  the 
infinite  integral  occurring  in  Eq.  (77) .   The  first  term  in 
the  integrand  is  singular  at  y  =  a .   Since  the  numerator 
approaches  zero  like  (y-a)  near  a,  Garrison  [Ref.  2]  sub- 
tracted l/(u-a)  from  the  singular  term  in  the  range  0  <_  y  <_  2a 
and  added  the  contribution  of  the  integral  of  1/ (y-a) (which 
in  this  case  was  zero  )  .  This  converted  the  singular  integrand 
to  a  regular  function.   Applying  this  technique  to  numeri- 
cally evaluate  the  resulting  integral  between  0  and  2a, 
Simpson's  three-eights  rule  is  used  with  an  odd  number  of 
subdivisions,  thus  ensuring  that  the  point  y  =  a  is  not 
encountered. 

With  f,  determined  from  the  inversion  of  Eq.  (99) ,  the 
first-order  potential  and  its  derivatives  may  be  evaluated 
on  the  cylinder  surface.   Replacing  the  surface  integrals 
with  summations,  Eq.  (85)  and  its  derivatives  may  be  written 
as 


Uli   =    6ij(v)flj  i,j    =    1'2'3"'"m  (100) 


uLi  =  wv)fij  (101) 


u?    •    =    S    •  •  (v)f.  •  (102) 

lyi        pyij  lj 
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s    s 
where  u. •  ,  u,  •,  etc.  denote  functions  evaluated  at  the 

ith  nodal  point  on  the  cylindrical  surface.   The  complex 

matrices  in  Eqs .  (100-102)  are  defined  as  follows: 

S.,(v)  =  i.     /G(x.  ,y.  ;£,n;v)dS    i,j  =  1,2, ...,m  (103) 

ID       2tt  ^  /    i   i         1 

Wv)  =^  q  ^x(x1.yii5.niv)dS1  (104) 

A^lj 

WV)  =  J?  Aq  yQ7(*i.y1l5.T„v)dS1  (105) 

Ablj 


The  first-order   incident  potential  and  its  derivatives  may 
be  evaluated  directly  at  each  nodal  point  after  differentia- 
tion of  Eq.  (30)  as  needed,  and  combined  with  the  scattering 
potential  and  its  respective  derivatives. 

To  solve  the  second-order  problem,  the  function  f*(0 
required  in  Eq.  (95), and  defined  by  Eq .  (52),  must  be 
evaluated  numerically.   Dividing  the  mean  free  surface 
(y  =  0)  in  the  vicinity  of  the  cylinder  into  n  equal  incre- 
ments, f*  is  evaluated  at  each  nodal  point.   The  numerical 
solution  again  uses  Eqs.  (100-105)  with  y.  =  0  for  the 
purpose  of  evaluating  u, •  and  its  derivatives  at  the  mean 
free  surface  nodal  points.   The  first-order  wave  potential 
and  its  derivatives  are  evaluated  directly  at  each  nodal 
point,  and  combined  with  the  scatter  potential  and  its 
derivatives  to  solve  Eq.  (52)  for  f*. 
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Having  determined  f*  at  all  nodal  points  on  the  free 

S* 

surface,  the  boundary-value  problem  for  u,,   given  by  Eqs. 

S*       S* 
(54-56)  may  be  solved,  i.e.  u0   and  un   may  be  evaluated 

2        2n 

on  the  cylinder  by  solution  of  Eqs.  (92)  and  (95)  respec- 
tively.  Expressing  the  integrals  as  numerical  summations  in 
complex  matrix  form  yields: 


u2ni  =  fiaij<4v>         i   =  1,2, ...,m  (106) 

j  =  1,2, . .  .  ,n 


u,*  =  fJBJ. (4v)  (107) 

z  J.      j  x  J 


where 


a*-(4v)  =i     /G*(x. ,y. ;?,n;4v)  d£ 
ID        tt  Ac  J  11 


(108) 


AS2j 


6ij(4v)  =  2T  AO  /GJ(xi^i'^'4^  ^ 


(109) 


AS  r 


s*     s* 

Thus,  using  f * ,  both  u^   and  u2   may  be  directly  evaluated 
at  each  cylinder  surface  nodal  point. 

To  solve  the  boundary-value  problem  for  the  second  part 

S° 
of  the  second-order  scatter  potential,  u-  ,  as  specified 

by  Eqs.  (57-60),  Eq.  (88)  must  be  evaluated.   However,  to 

S° 

evaluate  Eq .  (8  8)  for  u-  ,  the  second-order  source  strength 

function,  f2/  must  be  evaluated  by  numerically  solving 
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the  integral  equation,  Eq.  (91) .   Replacing  the  integral 
equation  with  a  complex  matrix  summation: 


f2i  +  f2jaij(4v)  =  2ki       i/j  =  1'2""'m  (110) 


where  a.j(4v)  is  defined  by  Eq.  (96), 


f2j  =  f2(Xj,yj)  (in) 

and, 

k.  =  1 [n  (x. ,y.)sinh[2a(y.+h)]  (112) 

1    sinh  (ah)  L  y   x   1 

i2ax. 


+ 


1  -^"^^  s* 

inx(xi,yi)cosh[2a(yi+h) ]  e      -  ^2n(xi'yi) 


Once  a.  .  (4v)  is  evaluated,  Eq.  (110)  may  be  solved  to  obtain 
the  second-order  source  strength  function,  f2,  at  each  nodal 
point  on  the  cylinder  surface.   Stating  Eq.  (88)  in  summation 
form, 


S° 
u2i  =  f2j6ij(4v)        i,j  =  1»2,3,...,m  (113) 


where  S..(4v)  is  defined  by  Eq.  (103).   Using  Eq.  (113), 

s°    1] 

Up   may  be  determined  at  each  nodal  point  on  the  cylinder 

S* 
surface,  then  combined  with  u?   and  the  second-order  inci- 
dent potential,  u_ ,  to  evaluate  the  total  second-order 


potential,  u2 . 


48 


Since  the  first-order  potential,  u. /  and  its  derivatives, 
u,   and  u,  ,  as  well  as  the  second-order  potential,  u^,    are 
now  known  at  each  nodal  point  on  the  surface,  the  first- 
order,  second-order  and  steady-state  force  coefficients  and 
phase  shift  angles  may  be  determined  using  Eqs .  (74-76).   As 
shown  earlier,  the  integrals  may  be  replaced  by  summations, 
using  nodal  point  values  and  the  arc  length  increment,  A0. 
Once  each  integral  is  evaluated,  then  the  force  coefficient 
becomes  the  absolute  value  or  modulus  of  the  complex  integral 
result  and  the  phase  shift  angle  is  the  angle  whose  tangent 
is  the  imaginary  part  over  the  real  part  of  the  complex 
integral  result. 

D.   COMPUTER  SOLUTION 

In  order  to  utilize  the  method  of  solution  described 
in  Sections  B  and  C,  a  computer  program  was  developed  to 
carry  out  the  indicated  calculations.   Although  accuracy  was 
a  prime  consideration,  an  equally  important  requirement  was 
to  minimize  computer  run  time.   To  achieve  this,  every  attempt 
was  made  to  utilize  symmetry  and  to  generate  certain  constants 
and  matrices  to  eliminate  redundant  computations. 

The  program  consists  of  a  main  program  that  flows  in  the 
order  of  computation  presented  in  Sections  B  and  C,  along 
with  four  subroutines.   Two  subroutines  to  solve  the  first- 
order  Green's  function,  using  both  forms  of  the  function, 
Eqs.  (77)  and  (82),  were  developed  by  Garrison  [Ref.  2]. 
These  subroutines,  GREEN  and  GREENS  respectively,  have  been 
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modified  to  include  calculation  of  the  first-order  scatter 
potential  derivatives,  evaluation  of  the  first-order  scatter 
potential  and  its  derivatives  on  the  mean  free  surface,  y  =  0, 

and  evaluation  of  the  modified  Green's  function,  G* ,  for 

S*      S  * 

the  determination  of  u~   and  u0   at  each  cylinder  nodal 

2        2n 

point. 

For  elements  of  the  a  and  6  matrices  corresponding  to 
small  values  of  the  parameter  (x-£) ,  GREEN  is  used  while  for 
larger  values  of  (x-£) ,  GREENS  is  used.   With  the  exception 
of  the  diagonal  elements  in  Eqs .  (96)  and  (103-105)  and  a  few 
surface  nodal  points  in  Eqs.  (108-109),  the  majority  of 
elements  of  a  and  3  are  calculated  by  GREENS.   This  is  most 
fortunate  as  the  series  form  converges  rapidly,  requiring  much 
less  computer  time  than  the  equivalent  integral  form. 

Subroutine  GEODAT  reads  the  input  geometrical  data, 
generates  the  matrices  h.  and  k.  as  defined  by  Eqs.  (97)  and 
(112)  respectively,  and  calculates  certain  geometrical 
parameters  and  matrices  for  repeated  use  in  GREEN  and  GREENS. 
Subroutine  COMAT  inverts  complex  matrix  equations  and,  there- 
fore is  used  to  invert  Eqs.  (99)  and  (110)  to  determine  f,  and 
fy   respectively. 

A  cross-reference  between  text  and  computer  program 
nomenclature  is  given  in  Table  I. 
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TABLE  I:   Computer  Program  —  Text  Symbol  Cross-Ref erence 


Text 


Computer 
Program 


Text 


Computer 
Program 


a 

A 

d 

D 

fl 

F(I,1) 

f2 

F1(I,1) 

f* 

FS(L) 

Fll 

Cl(l) 

F12 

CI  (2) 

F21 

C2(l) 

F 
22 

C2(2) 

Fss 

21 

C3(l) 

Fss 

22 

C3(2) 

G 

GIJ,GIJEXT 

G* 

GIJ,GIJEXT 

G 
n 

GNIJ,GNJI 

G 

X 

GXIJ,GXJI 

Gy 

GYIJ,GYJI 

G 

yy 

GYY 

h 

H 

2hi 

HH(I,1) 

ki 

PK(I,1) 

m 

NPTS 

n 

NSPTS 

n 


n 


u- 


u 


u 


lx 


ly 


u..  (surface) 


u,   (surface) 


ulv  (surface) 


ANX(I) 
ANY  (I) 
U1(I) 
U1X(I) 
U1Y(I) 
U1IS 

U1ISX 

U1ISY 


ulyy  (sur 

face) 

U1ISYY 

u,  (surf 

ace) 

U1SS 

u,   (sur 
lx 

face) 

U1SSX 

uly  (sur 

face) 

U1SSY 

ulyy  (sur 

face) 

U1SSYY 

u2 

U2(I) 

X 

X(I) 

y 

Y(I) 

a 

ALPHA (I, J) 

B 

BETA (I, J) 

ex 

BETAX(I,J) 

6y 

BETAY(I,J) 

A6  (cylin 

der) 

DELTHE 

A£  (surface) 

DELX 
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Text 


Computer 
Program 


Text 


Computer 
Program 


11 

PHASEl(l) 

V 

12 

PHASEK2) 

yk(v) 

21 

PHASE2(1) 

yk(4v) 

22 

PHASE2(2) 

V 

5 

X(J) 

4v 

n 

Y(J) 

■n 

AMU 
AMU(K) 
AMU4 (K) 

ANU 

ANU4 
PI 


c;? 


IV.   DISCUSSION  AND  RESULTS 

A.   SELECTION  OF  COMPUTER  PROGRAM  PARAMETERS 

The  computer  program  was  developed  on  the  IBM-360 
computer  using  FORTRAN  H.   All  subroutines  were  compiled  on 
DATA  CELL  to  minimize  compile  time  and  core  size.   In  order 
to  maximize  accuracy  while  minimizing  computational  time, 
two  key  computer  parameters  were  evaluated  to  determine 
optimum  values:  cylinder  nodal  points  and  free  surface  model 
points. 

Using  the  first-order  solution  only,  the  number  of  nodal 
points  on  the  cylinder  surface  was  varied  from  12  to  60  and 
compared  with  the  results  determined  by  Garrison  [Ref.  2]. 
The  minimum  number  of  cylinder  surface  nodal  points  that 
provided  accuracy  within  one  percent  of  those  obtained  by 
Garrison  using  60  points,  was  24.   Accordingly,  the  good 
correlation  from  60  points  down  to  24  points  led  to  the 
selection  of  24  nodal  points  on  the  cylinder  surface. 

Since  the  free  surface  integral  was  to  be  carried  out 
from  -«  to  +°°,  the  next  step  was  to  establish  the  finite 
interval  of  convergence  on  the  surface  as  well  as  the  sub- 
division size.   After  the  outer  limits  of  the  free  surface 
integral  were  determined,  the  subdivision  size  was  varied 
from  4  to  128  subdivisions  per  second-order  wave  length 
holding  the  total  interval  constant.   Above  64  subdivisions 
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the  results  did  not  vary  more  than  one  percent,  leading  to 
the  selection  of  64  subdivisions  per  second-order  wave  length 
on  the  free  surface  for  purposes  of  making  calculations. 

Additionally,  the  total  interval,  2 A,  was  varied  from 
one  to  eight  wave  lengths  in  both  the  positive  and  negative  x 
directions.   The  second-order  force  coefficients  and  their 
respective  phase  shift  angles  converged  to  a  small  value 
varying  periodically.   Convergence  to  this  limit  cycle 
occurred  between  the  first  and  second  wave  length  from 
the  center  of  the  cylinder.   Therefore,  a  total  interval  of 
from  -2  to  +2  second-order  wave  lengths  was  chosen  as  adequate 
for  generation  of  the  numerical  results. 

Although  the  values  of  the  second-order  force  coefficients 
and  their  respective  phase  shift  angles  tended  to  converge 
with  increasing  the  limits  of  the  free  surface  integration, 
there  remained,  in  general,  the  small  limit  cycle  as  noted 
above.   The  amplitude  of  the  cycle  was  a  function  of  the 
parameter,   a,  and  was  of  significant  magnitude  only  for 
values  of  a  less  than  0.25.   The  effects  of  a  on  the  limit 
cycle  for  one  of  the  second-order  parameters,  horizontal  force 
coefficient,  is  demonstrated  in  Fig.  (2).   This  represents  a 
plot  of  percent  variation  of  the  force  coefficient  from  the 
mean  versus  the  surface  interval  size,  X/(L/2),  corresponding 
to  values  of  a  of  0.25,  0.20  and  0.15.   Figure  (2)  is 
representative  of  the  behavior  of  all  second-order  force 
coefficients  and  phase  shift  angles  and  demonstrates  the 
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increase  in  the  amplitude  of  the  variation  with  a  decrease 
in  a.   The  amplitude  of  the  periodic  variation  ranges  from 
up  to  fifty  percent  at  a  wave  number  of  0.15  and  depth  of 
d  =  1.4,  to  less  than  two  percent  for  all  values  of  wave 
number  above  0.25. 

B.   RANGE  OF  APPLICABILITY 

In  a  complex  computer  program  of  the  type  developed  to 
carry  out  the  present  numerical  scheme,  certain  limits  with 
respect  to  numerical  stability,  etc.  show  up.   These  computing 
limits  arise  for  various  reasons  such  as  overflows  caused  by 
computing  hyperbolic  functions  of  very  large  numbers, 
numerical  convergence  instabilities,  etc.   While  no  attempt 
was  made  to  investigate  each  of  these  numerical  problems ,  a 
list  of  the  limitations  of  the  particular  program  developed 
in  this  work  to  carry  out  the  computations  are  as  follows: 

minimum  a  —  dependent  upon  acceptable  error,  but  0.25 
or  less 

maximum  a  —  deep  water  condition  reached 

minimum  h  —  3 

maximum  h  —  deep  water  condition  reached  for  a  >  0.25, 
h  =  20  for  a  <  0.25 

minimum  d  —  from  1.4  at  h  equal  to  or  less  than  5  down 
to  2.5  at  h  =  20 

minimum  SMIN  —  0.12  when  cylinder  depth  is  other  than 

near  the  free  surface 

maximum  SMIN  —  0.30  when  cylinder  depth  is  near  the 

free  surface 
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C.   RESULTS 

A  representative  set  of  results  for  the  first-order  and 
second-order  force  coefficients  have  been  generated  for  a 
water  depth  of  h  =  5.0,  submergence  depths  of  d  =  1.5,  2.0, 
2.5,  and  3.0,  with  a  ranging  from  0.15  to  1.2.   All  of  the 
numerical  results  presented  herein  were  based  on  24  sub- 
divisions on  the  circular  cylinder  and  64  subdivisions  per 
second-order  wave  length  on  the  free  surface  integration. 
The  surface  integral  was  carried  out  from  x  =  -L  to  +L  since 
this  was  found  to  be  adequate  for  convergence.   Since  the 
second-order  wave  length  is  half  the  first-order  wave  length, 
L,  the  total  interval,  is  equal  to  four  second-order  wave 
lengths. 

The  first-order  horizontal  and  vertical  force  coefficients 
are  presented  in  Figs.  (3)  and  (4),  respectively.   It  may  be 
noted  that,  in  general,  the  forces  decrease  with  depth  of 
submergence  according  to  expectation  since  the  wave  action 
dies  out  with  depth. 

The  second-order  horizontal  and  vertical  force  coefficients 
are  presented  in  Figs.  (5)  and  (6),  respectively.   Generally, 
the  results  show  the  second-order  effect  to  be  relatively  more 
important  as  the  cylinder  approaches  the  free  surface.   It 
might  be  suspected  from  this  that  the  second-order  contribution 
would  be  much  more  significant  in  the  case  of  surface  piercing 
or  floating  bodies . 
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The  horizontal  and  vertical  steady-state  force  coefficients 
are  shown  in  Figs.  (7)  and  (8) ,  respectively.   These  results 
show  that  the  horizontal  steady-state  force  coefficient  is  very 
small  in  general.   The  horizontal  force  can  be  shown  to  be 
proportional  to  the  momentum  flux  of  the  reflected  wave  and 
the  reflected  wave  is  small.   In  fact,  in  the  case  of  infinite 
water  depth,  Dean  [Ref.  1]  showed  that  the  reflection  was 
exactly  zero,  and,  accordingly,  the  steady-state  horizontal  force 
is  likewise  zero. 

It  may  be  noted  that  the  steady-state  vertical  force  is 
positive.   This  results  from  the  fact  that  the  velocities  are 
largest  on  the  top  of  the  cylinder  and,  accordingly,  the  average 
pressure  is  reduced  on  the  top  of  the  cylinder. 

The  phase  shift  angles  of  the  first-order  and  second-order 
forces  are  shown  in  Figs.  (9)  -  (12) . 

To  demonstrate  the  effect  that  the  second-order  terms  have 
on  the  forces  on  the  cylinder  and  their  respective  phase  shift 
angles,  a  comparison  was  made  between  the  first-  and  second- 
order  results.   Figures  (13)  -  (13)  are  plots  of  the  horizontal 
and  vertical  dimensionless  forces  on  the  cylinder  versus  time 
over  one  complete  wave  cycle  for  both  the  first-order  solution 
alone  and  the  combined  first-order  and  second-order  effects. 
Additionally,  a  plot  of  the  incident  wave  is  included  to 
demonstrate  the  phase  shift  involved  (the  amplitude  of  the 
incident  wave  has  no  significance).   A  mean  water  depth,  h,  of 
5.0,  a  cylinder  depth,  d,  of  1.5  and  a  wave  height,  H,  of 
0.5  were  used,  with  the  wave  number,  a,  assigned  three  values, 
0.25,  0.5  and  1.0. 
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V.   CONCLUSIONS 

A  computer  program  has  been  developed  to  carry  out  the 
numerical  solution  of  the  second-order  wave  —  cylinder 
interaction  problem.   The  first-order,  second-order,  and 
steady-state  force  coefficients  were  determined  for  the 
submerged  horizontal  cylinder. 

An  approximate  method  for  dealing  with  the  second-order, 
nonhomogeneous  free  surface  boundary  condition  was  developed 
which  appears  to  converge  except  at  small  values  of  a,  i.e., 
at  very  large  wave  lengths. 
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APPENDIX  A 

COMPUTER  PROGRAM  LISTING 

THIS  PROGRAM  CALCULATES  WAVE  FORCES  AND  THEIR  PHASE 

SHIFT  ANGLES  FOR  WAVE  INTERACTION  WITH  A  SUBMERGED 

HORIZONTAL  CYLINDER 

A  =  2*PI*CYLINDER  RADIUS/WAVE  LENGTH  =  SIGMA**2* 

CYLINDER  RADIUS/ACCELERATION  OF  GRAVITY 

H  =  MEAN  WATER  DEPTH/CYLINDER  RADIUS 

D  =  CYLINDER  DEPTH/CYLINDER  RADIUS 

NPTS  =  NUMBER  OF  CYLINDER  SURFACE  ELEMENTS  AND  NODAL 

POINTS  FOR  NUMERICAL  EVALUATION 

NSPTS  =  NUMBER  OF  FREE  SURFACE  ELEMENTS  AND  NODAL 

POINTS  FOR  NUMERICAL  EVALUATION 

COMPLEX  GIJ,GNI J,GNJI ,GXI J, GX J  I , GYI J, GY J  I , G IJ EXT , GYY, 
1DET,SUM1,SUM2,SUM3,SUM4,U1I  S,U1ISX,U1ISY,U1ISYY, 
2U1SS,U1SSX,U1SSY,U1SSYY 

COMPLEX  ALPHA  (24,  24)  , BET A (24, 24)  , BET AX ( 24, 24)  , 
1BETAY(2  4,24) ,HH( 24,1) ,F(24, 1) ,PK(24, 1) ,F1(24,1) , 
2FS(500) ,U1(24) ,U1X(24) ,U1Y( 24) ,U2(24) ,J2SC1(24) , 
3U2SCNK24)  ,U2SCO(24) 

DIMENSION  X<24) ,Y(24)  ,ANX( 24) , ANY ( 24) , CHY I  2 5 ) ,CHY2( 25) 
1,SHY(25) ,SHY2(25)  ,COEFG(200) ,COEFG2(200)  , 
2COSAMU(  200,2  5) ,C0SAM2( 200,25) , SINAMU(2J0  ,25), 
3SINAM2( 200,25) , AMU (200) ,AMU4(200) ,SH2Y( 24) ,CH2Y( 24) 

DIMENSION  CI (2) ,C2(2) ,C3(2) ,PHASE1(2) ,PHASE2( 2) 

COMMON/CPX/HH,PK 

CCMMON/VAR/X,Y,ANX,ANY 

C0MMCN/GSHY/ChY,CHY2,SHY,SHY2,SH2Y,CH2Y 

CCMM0N/GMU/C0SAMU,C0SAM2,SINAMU,SINAM2, AMU , AMU4 , COEFG, 
1C0EFG2 

COMMON/CONST/H,D,DELTHE,SMI N, NPTS, NSPTS 

COMMON/CPI/PI 

EQUIVALENCE  (Hh(  1),F(  1)  ) 

EQUIVALENCE  ( PK (  1  )  ,  F  1  (  1 )  ) 

PI=3. 141592 

READ    INPUT    DATA    AND    CALCULATE    REQUIRED    REPEATING 
DATA    AND    ARRAYS 

CALL    GEODAT    ( A, A2 , ANU, ANU4, SH2 AH , SH2AH2 , SHAH, SHAH2 , 
1CHAH,CHAH2, AO, A A , 3B , CC , DD , A02 , A A2 , 8B2 , CC 2, 0D2 ) 

DO    100     1=1, NPTS 

DO    100    J=1,I 

XV=X( I ) 

YV=Y( I) 

XC=X( J) 

YC=Y( J) 

ANXI=ANX( I ) 

ANYI=ANY( I ) 

ANXJ=ANX( J) 

ANYJ=ANY( J ) 

IF(ABS(X( I)-X( J) ) .LT.SMIN)     GO    TO    50 

SHYI=SHY( I ) 

CHYI=CHY( I ) 

SHYJ=SHY( J) 

CHYJ=CHY( 


CALL  GREE 
1C0EFG,SHY 
2ANXJ, ANYJ 

GC    TO    90 

50    CONTINUE 


(J) 

ENS     I  A, ANU, S H2 AH, SH AH, C HAH, COS AMU, SIN AMU, AMU, 
YI ,CHYI,SHYJ,CHYJ, I,J,XV,YV,XC,YC,ANXI,ANYI, 
J,GIJ,GNIJ,GNJI,GXIJ,GXJI,GYIJ,GYJI ,GYY) 


CALL    GREEN    ( A, ANU , SH2 AH, SHAH,CHAH, AO, AA , BB , CC , DD , I , J , 
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1XV,YV,XC,YC,ANXI,ANYI ,GI J , GNI J , GXI J, GYI J ,GYY) 

IF(I.EQ.J)     GNJI=GNIJ 
IF    (I.EQ.J)     GXJI=GXIJ 
IF    (I.EQ.J)     GYJI=GYIJ 
IF( I.EQ.J)     GO    TO    90 

CALL    GREEN    ( A , ANU, SH2 AH, SHAH,ChAH, AOt AA , BB , CC ,DD t J » I » 
lXCfYCi XVfYVi ANXJiANYJ,GI JEXTtGNJI.GXJIf GYJI ,GYY) 

90    CONTINUE 

CALCULATE  THE  FIRST-ORDER  ALPHA  AND  BETA  MATRICES 

ALPHA (  I,J)=(  l./PI )*GNIJ*DELTHE 

ALPHA! J, Ii=( l./PI )*GNJI*DELTHE 

BETA( I,J)=(1./(2.*PI) )*GI J*DELTHE 

BETA( J,  U=BETA(I  ,  J) 

3ETAX(I,J)=(1./(2.*PI) )*GXI J*DELTHE 

BETAX( J, I)=( l./( 2.*PI ) )*SXJI*DELTHE 

BETAY( I,J)=(1./(2.*PI) )*GYI J*DELTHE 

BETAY(J,I)=(1./(2.*PI) )*GYJI*DELTHE 
100  CONTINUE 

DO  120  I=1,NPTS 

ALPHA { 1,1) = ALPHA* I , I )+CMPLX( 1.0,0.0) 

BETAX( I  ,  I  )=BETAX( I  ,  I) *ANX( I ) *CMPLX (0 .5 , 0 .0) 

BETAYU  ,I)=BETAY(  I  , I)+ANY(  I  ) *CMPLX (0 .5 , 0 . 0) 
120  CONTINUE 

GENERATION  OF  THE  FIRST-ORDER  DISTRIBUTION  FUNCTION, 
F(I,1),  BY  INVERSION  OF  ALPHA{  I , J  )  *F( I , 1 )  =  HH(I,1) 

CALL  COMAT  (  24,  1  ,  ALPHA  ,HH  ,DE  T,  IND  ICA  ) 

THE  HH(I,1)  MATRIX  IS  REPLACED  BY  THE  DISTRIBUTION 
FUNCTION, F( I ,1 ) 

COMPUTATION  OF  THE  FIRST-ORDER  SCATTERING  POTENTIAL 
FUNCTION,  U1SC(I),  AND  ITS  X  AND  Y  PARTIAL  DERIVATIVES 
UISCX(I)  AND  U1SCYII),  AND  COMBINATION  WITH  THE 
INCIDENT  POTENTIAL  COUNTERPARTS  TO  COMPUTE  THE  TOTAL 
POTENTIAL  FUNCTION  AND  ITS  X  AND  Y  DERIVATIVES 

DO  155  1=1,NPTS 

SUM1=(0. 0,0.0) 

SUM2=(0. 0,0.0) 

SUM3=(0. 0,0.0) 

DO  150  J=1,NPTS 

SUM1=SUMH-F(  J,l)  *BETA(  I,J) 

SUM2=SUM2+F( J,1)*BETAX( I, J) 

SUM3=SUM3+F( J,1)*BETAY( I, J) 
150  CONTINUE 

U1(I)=SUM1-(CHY( I )/( A*CHAH) ) *CEX P{ CMPLX < 0 . 0 , A*X ( I ) ) ) 

U1X(I)=SUM2-CMPLX(0.0,1.0)*CHY( I ) *CEXP ( CMPLXt 0. 0 , 
1A*X( I ) ) J/CHAH 

U1Y(I)=SUM3-SHY(  I  )  *C  EX  P(  CMPLX(  0  .  0  ,  A*X  (  I  )  )  J/CHAH 
155  CONTINUE 

EVALUATION  OF  THE  FREE  SURFACE  PRESSURE  DISTRIBUTION 
SOURCE  STRENGTH  FUNCTION,  FS(L) 

DELX=0.015625*PI/A 
WRITE  (6,180)  DELX 
180  FORMAT  ( 5X//5X , • DELX=« , Fl 0. 5/// ) 
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SP=0.0 

CST=.5*DELTHE/PI 

IC0UNT=1 

DO    250    L=1,NSPTS 

SUM1=(0.0»0.0) 

SUM2=(0. 0,0.0) 

SUM3=(0. 0,0.0) 

SUM4=(0. 0,0.0) 

DO    245     I=1,NPTS 

XV=SP 

YV=0.0 

XC=X(I) 

YC=Y(I ) 

ANXI=0.0 

ANYI=1.0 

ANXJ=ANX(  I  ) 

A.NYJ=ANY(  I  ) 

SHYI=SHY(25) 

CHYI=CHY(25) 

SHYJ=SHY( I) 

CHYJ=CHY( I ) 

IF     (ABS(XV-XC) .LE.SMIN)    GO    TO    240 

CALL  GREENS  ( A , ANU , SH2 AH, SHAH, CHAH,COSAHU, S INAMU , AMU, 
1C0EFG,SHYI ,CHYI, SHY  J , CHY J , 25 , I ,XV , YV ,XC , YC , ANXI , ANYI , 
2ANXJ,ANYJ,GIJ,GNIJ,GNJI,GXI  J  ,  GX  J  I  ,GYIJ  ,  GY  JI  ,GYY) 

GO    TO    241 

240  CONTINUE 

CALL    GREEN    ( A , ANU, SH? AH , SHAH,CH AH, AO, AA , BB , CC , DO, 2 5, I , 
1XV,YV,XC,YC,ANXI , ANY  I ,GIJ , GNI J , GX I  J, GY IJ , GY Y) 

241  CONTINUE 

SUM1  =  SUMH-CST*GIJ*F(  1,1) 

SUM2-SUM2+CST*SXIJ*F(  1,1) 

SUM3=SUM3+CST*GYIJ*F( 1,1) 

SUM4=SUM4+CST*GYY*F( 1,1) 
245  CONTINUE 

U1IS=(-1./A)*CEXP(CMPLXO.O,A*SP)  ) 

U1ISX=CMPLX(0.0,-1.0)*CEXP(CMPLX(0.0,A*SP) ) 

U1ISY=-TANH<A*H)*CEXP(CMPLX(0.0,A*SP) ) 

U1ISYY=-A*CEXP(CMPLX(0.0,A*SP)  ) 

U1SS=SUM1 

U1SSX=SUM2 

U1SSY=SUM3 

U1SSYY=SUM4 

FS(L)  =  ( 2.*A/{3.*ANU) ) * ( Ul I S*U1 SSYY HJ1 I S YY*U1SS+U1  SS* 
1U1SSYY-6.*U1ISY*U1SSY-3.*U1SSY*U1SSY-4.*UIISX*U1SSX 
2-2.*UlSSX*UlSSX) 

IF  ( ICOUNT.EQ.O)  GO  TO  248 

SP=FL0AT(L*l)*DELX/2. 

ICOUNT=0 

GO  TO  250 
248  SP=-SP 

IC0UNT=1 
250  CONTINUE 

CALCULATION  OF  THE  PARTICULAR  SOLUTION  PORTION  OF  THE 
SECOND-ORDER  SCATTERING  POTENTIAL  AND  ITS  NORMAL 
DERIVATIVE,  U2SCKI)  AND  U2SCNKIJ 

CST=.5*DELX/PI 
DO  350  I=1,NPTS 
SUM1=(0. 0,0.0) 
SUM2=(0. 0,0.0) 
SP=0.0 
IC0UNT=1 

DO  345  L=1,MSPTS 
XV=X( I ) 


YV  =  YU  ) 
XS=SP 
YS=0.0 
ANXI=ANX( I) 
ANYI=ANY(I ) 
ANXJ=0.0 
ANYJ=1.0 
SHYI=SHY(I ) 
CHYI=CHY(I ) 
SHYJ=SHY(25) 
CHYJ=CHY(25) 
IF    (ABS(XV-XS) 


LE.SMIN)    GO    TO    340 


CALL    GREENS     ( A2 , ANU4 , SH2AH2 , SHAH2 ,CHAH2 , C0SAH2, SI NAM2, 
1AMU4,C0EFG2,SHYI,CHYI,SHYJ,CHYJ, I,25,XV,YV,XS,YS, 
2ANXI,ANYI,  ANXJ,ANYJ,GU,GNI  J,GNJI,GXI  J,GXJI ,GYU,GYJI, 
3GYY) 


GO  TO  341 
340  CONTINUE 


CALL  GREEN  ( A2 » ANU4, SH2 AH2 , SHAH2 , CHAH2, A02, AA2 , BB2 , 
Lt 002i I , 25,XV,YV,XS,YS,ANXI , ANY  I ,GU , GNI J , GX IJ , GY IJ , 
2GYY) 
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341    CONTINUE 

SUM1=SUM1+CST*GIJ*FS( L) 

SUM2=SUM2+CST*GNIJ*FS(L) 

IF    ( ICOUNT.EQ.O)     GO    TO    290 

SP=FL0AT(L+l)*DELX/2. 

ICOUNT=0 

GO    TO    345 
290    SP=-SP 

IC0UNT=1 
345    CONTINUE 

U2SCK  I  )=SUM1 

U2SCNK  I)=SUM2 
350  CONTINUE 

CALCULATION  OF  THE  HOMOGENEOUS  SOLUTION!  PORTION 
OF  THE  SECOND-ORDER  SCATTERING  POTENTIAL,  U2SC0U) 

DO  500  I=1»NPTS 
DG  500  J=1,I 
XV=X( I ) 
YV=Yl I ) 
XC=X( J) 
YC=Y( J) 
ANXI=ANX( I ) 
ANYI=ANY( I ) 
ANXJ=ANX( J) 
ANYJ=ANY( J ) 
IF(ABS(X( I)-X( J)  ) 
SHYI=SHY(I ) 
CHYI=CHY( I ) 
SHYJ=SHY{ J) 
CHYJ=CHY( J ) 

CALL  GREENS  ( A2 , ANU4, SH2AH2 , SHAH2 ,CHAH2 , COS AM2» S I NAM2t 
1AMU4,C0EFG2»SHYI T CHY I t SHY  J , CHY J , I , J, XV, Y V ,XC, YC , ANX I , 
2  ANY  I  ,  ANX  J  ,  ANY  J,  G  I  J,  GNU,  GN  J  I ,  GX  IJ  ,  GX  J  I ,  GYIJ  ,GYJ  I  ,GYY) 

GC  TO  490 
450  CONTINUE 

CALL  GREEN  ( A2 , ANU4, SH2AH2 , SHAH2 ,CHAH2 , A02 , AA 2 ,BB2 ,CC2 
1  ,002  t  I  ,  J,XV,YV,XC,YC,ANXI  ,ANYI  ,  GIJ  ,  GNU  ,  GX  IJ,  GY  I  J  , 
2GYY) 


LT.SMIN)  GO  TO  450 


IF(I.EQ.J)  GNJI=GNU 
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IF< I.EQ.J)  GO  TO  490 

CALL  GREEN  ( A2 , ANU4, SH2AH2 , SHAH2 , CHAH2 , A02 , AA 2 , BB 2 TCC2 
It  DD2i  Jt  IfXCYCXVtYVfANXJt  ANYJ,GIJEXT,GNJIi  GXJI»GYJI, 
2GYY) 

490  CONTINUE 

ALPHA (I ,  J)=(1./PI )*GYIJ*DELTHE 

ALPHA ( J,I)=< l./PI)*GYJI*DELTHE 

BETA( I,J)=(1./(2.*PI) )*GIJ*DELTHE 

BETA( J, I)=BETA( I , J) 
500  CONTINUE 

DO  510  I=1,NPTS 

PK( I,1)=2*(PK(I,1)-U2SCN1U )) 
510  CONTINUE 

DO  520  I=1,NPTS 

ALPHA (I , I)=ALPHA( I, I)+CMPLX( 1.0,0.0) 
520  CONTINUE 

CALL  COMAT  ( 24, 1 , ALPHA , PK , DET, IND ICA ) 

DO  540  I=1,NPTS 

SUM=(0. 0,0.0) 

DO  530  J=1,NPTS 

SUM=SUM+F1( J,1)*BETA{ I, J) 
530  CONTINUE 

U2SC0( I >=SUM 
540  CONTINUE 

EVALUATION  OF  THE  TOTAL  SECOND-ORDER  POTENTIAL,  U2(I) 

DO  550  1=1 ,NPTS 

U2(I)=U2SC1( I)+U2SC0( I)-(CHY(I ) /(2.*A*( SHAH**4) ) )* 
1CEXP(CMPLX(0.0,2.*A*X( I ) ) ) 
550  CONTINUE 

WRITE  (6,560) 
560  FORMAT  ( 5X//9X , • I «  ,  15X , • Ul ( I )  ' , 24X , • U1X (  I ) ■  ,25X , 
1«U1Y(I ) • ,24X, 'U2( I )'//) 

DO  580  I=1,NPTS 

WRITE  (6,570)  I , U  1  ( I ) , U1X ( I ) ,U1 Y ( I ) , U2( I ) 
570  FORMAT  (  5X  ,  I  5 , 4(  F 16 .6 , F  14.6  )  ) 
580  CONTINUE 

EVALUATION  OF  THE  FIRST-ORDER,  SECOND-ORDER  PERIODIC, 
AND  SECOND-ORDER  STEADY  STATE  FORCE  COEFFICIENTS  AND 
THE  PERIODIC  FORCE  PHASE  SHIFT  ANGLES 

SUM1=(0. 0,0.0) 

SUM2=(0. 0,0.0) 

SUM3=(0. 0,0.0) 

SUM4=(0. 0,0.0) 

SUM5=0.0 

SUM6=0.Q 

DO  720  I=1,NPTS 

SUM1=SUM1*U1 ( I)*ANX( I )*DELTHE 

SUM2=SUM2+U1( I}*ANY(I )*DELTHE 

SUM3  =  SUM3  +  ( (6.*ANU*ANU*U2(  I ) /A ) -U1X( I ) *U1 X (  I)-U1Y(  I  ) 
1*U1Y(I) )*ANX( I)*DELTHE 

SUM4=SUM4+( (6.*ANU*ANU*U2(  I ) /A )-UlX( I ) *U IX <  I)-U1Y( I ) 
1*U1Y( I)  )*ANY( I)*DELTHE 

SUM5  =  SUM5  +  (CABS(U1X(  I)*U1X( I ) ) +CA8S(U1Y(  I)*U1Y(I)  ) 
L+ANU*ANU/A**2-1 . 0)*ANX( I)*QELTHE 

SUM6=SUM6*(CABS(U1X(  I)*U1X(  I ) ) +CABS(U1Y(  I)*U1Y( I )  ) 
1*-ANU*ANU/A**2-1  .  0  )  *ANY(  I  )  *D  EL  THE 
720  CONTINUE 

C1(1)=CABS(SUM1*A) 

C1(2)=CABS(SUM2*A) 

C2(1)=CABS(SUM3*A*A/(4.*ANU) ) 
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C2(2)=CABS(SUM4*A*A/(4.*ANU) ) 

C3(1)=SUM5*A*A/(4.*ANU) 

C3(2)=SUM6*A*A/(4.*ANU) 

PHASE 1( 1)=ATAN2( AIMAG( SUMi*A) , RE AL( SUM1* A ) ) 

PHASE1(2)=ATAN2(  A IMAG< SUM2* A ) , REAL ( SUM2* A )  ) 

PHASE  2 ( 1)=ATAN2( A IMAG ( SUM3*A*A/ ( 4.*ANU )  )  , RE AL ( SUM3*A*A 
1/(4. *ANU) ) ) 

PHASE2  1 2)=ATAN2( AIMAG(SUM4*A*A/(4.*ANU) ) ,RE AL( SUM4*A*A 
1/(4. *ANU) ) ) 

WRITE    (6,730)    C 1 ( 1 ) , C 1 ( 2) , C 2 ( 1 ) , C2 ( 2) t C 3 ( 1 )  ,C3( 2 ) , 
1  PHASE  1 (1)  ,PHASE1 ( 2  ) , PH ASE2 ( 1 ) , PH ASE2( 2 ) 
730    FORMAT     (  5X// /5X ,  • CI ( 1 )  =  » , F8 . 5, 5X , • CI ( 2 )  =  • , F8. 5// 5X , 
1,C2(1)=,,F8.5,5X,'C2(2)=«,F8.5//5X,,C3(1)=',F8.5,5X, 
2»C3(2)=«,F8.5//5X,,PHASE1(1)=',F8.5,5X,1 PHASE1(2)=», 
3F3.5//5X, • PHASE  2 ( 1 )  =  »  , F8. 5 , 5X,  « PHASE2(2)=»  , F8.5//) 

STOP 

END 

THIS    SUBROUTINE    READS    THE    INPUT    GEOMETRICAL    DATA    AND 
CALCULATES    THE    FIRST    200    ROOTS    OF    AMU*T AN( AMU*H)     + 
ANU    =    0    AND    AMJ4*TAM( AMU4*H)     +    ANU4    =    0.        IT    ALSO 
GENERATES     ARRAYS    OF    CERTAIN    FUNCTIONS    AND    COEFFICIENTS 
USED    IN    GREEN    AND    GREENS    AS    WELL    AS    THE    HH    AND    PK 
MATRICES 

SUBROUTINE    GEODAT     ( A , A2 , ANU , ANU4, SH2AH , SH2AH2 , SHAH , 
1SHAH2,CHAH,CHAH2,A0,AA,B6,CC,DD,A02,AA2 , BB2 ,CC2, DD2 ) 

COMPLEX    HH(24,1) ,PK(24,1) 

DIMENSION    X(24) , Y I  24 ) , ANX ( 24 ) , ANY( 24) ,CHY(2  5)  ,CHY2(25) 
1, SHY (25) , SHY2( 25) , CCEFG< 2  00 ) , C0EFG2 ( 200 ) , AMU (200) , 
2AMU4(20  0) ,C0SAMJ(20C,25) ,C0SAM2( 200,25) , SINAMU( 200,25) 
3,SINAM2(2  0C,25) , SH2Y ( 24) , CH2Y ( 24)  ,XX( 24) 

CGMMON/CPX/HHfPK 

CCMMON/VAR/X,Y,ANX,ANY 

CCMM0N/GSHY/CHY,CHY2, SHY , SHY2, SH2Y , CH2Y 

C0MM0N/C-MU/C0SAMUfC0SAM2tSINAMUf  SINAM2,  AMU  ,  AMU4,  COEFG, 
1C0EFG2 

C0MMON/C0NST/H,D,DELTHE,SMIN,NPTS,NSPTS 

COMMON/CPI/PI 

RE AD (5, 5)     A,H,D, SM  IN, NPTS , NSPTS 
5    FORMAT     (4F10.5,2I4) 

ANU=A*TANH(A*H) 

SH2AH=S1NH(2.*A*H) 

CHAH=C0SH( A*H) 

SHAH=SINH( A*H) 

AO=TANH(A*H) 

AA*1./CHAH**2 

BB=-H*SHAH/ICHAH**2) 

CC=-H*H*( l.-SHAri**2)/(3.*CHAH**3) 

D0=H**3*SHAH*( 2.*CHAH**2+3.*( 1.- SHAH** 2) ) / < CHAH**4* 12) 

DELTHE=2.*PI/NPTS 

THETA=0ELTHE/2. 

DO    15    1=1, NPTS 

X(I)=COS(THETA) 

Y(  I)  =  SIHTHETA)-D 

ANX(I)=X( I ) 

ANY( I)  =  Y( I  )*D 

SHY( I )=SINH( A*(H+Y( I ) ) ) 

CHY(I )=COSh( A*(H+Y( I ) ) ) 

SH2Y( I)=SINH(2.*A*(H+Y( I ) ) ) 

CH2Y( I )=C0SH(2.*A*(H+Y( I) ) ) 

HH( I,1)  =  (2./CHAH)*CMPLX(ANY(  I ) *SHY (  I )  , ANX (  I  )*CHY<  I )  )* 

1CEXP(CMPLX(0.0,     A*X(I)))  ,  „.  ,     »  » 

PK( I,1)=(ANY(  I)*SH2Y(I  )  +CMPLX ( 0  . 0  , 1 . 0 ) * ANX ( I)*CH2Y( I )) 
1*CEXP(CMPLX(0.0,2.0*A*X(I )) ) /SHAH**4 

THETA=THETA+DELTHE 
15    CCNTINUE 

SHY(25) =SHAH 

CHY(25)=ChAH 
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B=ANU*H 
DO    50    K  =  l, 200 
XX( 1)=PI*K 
DO    25    1=1,20 
1 1  =  1 

YY=XX( I ) 

XXU+1)=XX(1)-ATAN2{3,YY) 

IF(ABS( (XX( I+1)-XX( IJ )/XX( 1+1) ) .LT..0001)    GO    TO    30 
25    CONTINUE 

IF    ( II.GE.20)    GO    TO    27 
GO    TO    30 

27  WRITE    (6,28) 

28  FORMAT     {5X,'AMU    ROOT    DOES    NOT    CONVERGE'//) 
30    CONTINUE 

AMU(K)=XX(I I )/H 

C0EFG(K)=2.*PI*( AMU<K)*AMU< K) +ANU*ANU) /(  ANU*AMU( K) -H* 
1AMU<K)*(AMU(K)*AMU(K>+ANU*ANU) ) 

NPTS1=NPTS+1 

DO    40    I=1,NPTS1 

IF    (I  .EC.NPTS1 )     GO    TO    35 

COSAMU(K,I  )=COS(  AMU(K)*IH«-Y(  I)  )  ) 

SINAMU(Kt  I)  =  SIN(AMU(K)*<m-Y<I  )  )  ) 

GO    TO    4C 
35    COSAMU(K,I )=COS(AMU(K)*H) 

SINAMU(K,I )=SIN( AMU(K)*H) 
40    CONTINUE 
50    CONTINUE 

ANU4=4.*ANU 

B2=ANU4*H 

DO    70    K=l,200 

xxa  )  =  pi*k 

DO    55    1=1,20 
I  1  =  I 

YY=XX(  I  ) 

XX{H-1)=XXU)-ATAN2<B2,YY) 

IF(A8S( (XX( I+1)-XX( I ) )/XX( 1+11) .IT..0001 )     GO    TO    60 
55    CONTINUE 

IF    {  II.GE.20)     GO    TO    57 
GO    TO    6C 

57  WRITE    (6,5S) 

58  FORMAT     (5X,'4MU4    ROOT    DOES    NOT    CONVERGE'//) 
60    CONTINUE 

AMU4CK)  =XX(  I  I)/H 

C0EFG2(K)=2.*PI*( AMU4C K) *AMU4( K) +ANU4*ANU4) /( ANU4* 
1AMU4(K) -H*AMU4(K)*( AMU4 ( K) * AMU4 ( K) +ANU4* ANU4) ) 

NPTS1=NPTS+1 

DO    65    I=1,NPTS1 

IF    (I  .EQ.NPTS1 )    GO    TO    68 

C0SAM2(K,I )=COS( AMU4 ( K ) * ( H+Y ( I ) ) ) 

SINAM2(K, I )=SIN( AMU4( K )* ( H  +  Y( I ) ) ) 

GO  TO  65 
68  C0SAM2( K,I )=COS( AMU4(K)*H) 

SINAM2(  K,  I)=SIN( AMU4(K)*H) 
65  CONTINUE 
70  CONTINUE 

XX(1)=ANU4 

DO  80  1=1,20 

I  I  =  I 

Y2=XX( I ) 

XX(I+1)=4.*ANU/TANH(Y2*H) 

IF  (A3S( iXX(  I+1)-XX( I ))/XX(  1+1)  )  .LT. 0.0001)  GO  TO  85 
80  CONTINUE 

WRITE  (6,32) 
82  FORMAT  <5X,'A2  ROOT  DOES  NOT  CONVERGE'//) 
85  CONTINUE 

A2=XX(I I) 

SH2AH2=SINH(2.*A2*H) 

CHAH2=C0SH<A2*H) 

SHAH2=S INH( A2*H) 

A02=TANH( A2*H) 

AA2=1./CHAH2**2 

BB2=-H*SHAH2/(CHAH2**2) 
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CC2=-H*H*( l.-SHAH2**2)/(3.*CHAH2**3) 
DD2=(H**3)*SHAH2*( 2. *CHAH 2**2+3 . *( l.-SHAH2**2))/ 
1 ( (CHAH2**4)*12. ) 
DO    90    I = 1 , N P T S 
SHY2( I)=SINH(A2*(H+Y( I) )) 
CHY2(I )=COSH(A2*(H+Y( I ) ) ) 
90    CONTINUE 

WRITE    (6,95)     A,A2,D,H  ,ANU ,A NU4, NPTS, SMI N , NSPTS 
95    FORMAT     ( 5X///3X, ■ A  =  « , F10. 5 , 4X , ' A2  =  ' ,F ID . 5 ,4X, • D= • , 
1F10.5,4X,«H=,,F10.5,4X, •ANU=» , F 1 0 . 5 , 4X , ■ ANU4=« , F 1 0 . 5/ 
24X,«NPTS='  , I3,4X," SMIN=« ,F10.5,4X, 'NSPTS  =  » ,  13///) 
RETURN 
END 

THIS    SUBROUTINE    CALCULATES    G    AND    ITS    DERIVATIVES, 
GX,GY,GN,     AND    3YY,    3Y    USE    OF    THE     INTEGRAL    FORM    FOR    THE 
CASE    (X(I)    -    X(J))     LESS    THAN    SMIN 

SUBROUTINE    GREEN     ( A , ANU, SH2 AH, SH AH,CHAH , AO, AA ,88 , CC ,DD 
1,I,J,X,Y,XI , ETA, ANX,ANY,G,GN,GX,GY,GYY) 

COMPLEX  G,GN 

COMPLEX  GX,GY,GYY 

DIMENSION    TEST! 2  00) ,TESTT( 100) , TESTTT< 100) ,SUMOX( 15) , 
1SUM0Y(15) ,NNN(15) , TESTY Y( 200) 

CCMMON/CONST/H,D,DELTHE,SMIN,NPTS,NSPTS 

COMMON/CPI/PI 

P1(Y,ETA,XX,AMU)=C0SH(AMU*(Y+H) )*COSH(AMU*( ETA+H) )* 
1C0S( AMU*XX)/(COSH( AMU*H)**2 ) 

P2X(Y,ETA,X,XI,AMU)=-AMU*C0SH( AMU*(ETA+H) )*SIN( AMU* 
1 (X-XI ) )*COSH( AMU*(Y+H) ) / { COSH{ AMU*H ) **2 ) 

P2Y(Y,ETA,X,XI , AMU ) =AMU*COSH( AMU* (ETA+H) ) *S INH( AMU* 
1 (Y+H) ) *COS(AMU*(X-XI) ) / (COSH( AMU*H )**2 ) 

P2YY(Y,ETA,X.XI , AMU ) = AMU*AMU*COSH( AM J* ( Y +H ) )*CCSH( AMU* 
1( ETA+H) )*COS( AMU*( X-XI) ) / ( CCSH ( AMU*H ) **2 ) 

31(Y,ETA,XX,AMJ)=-P1(       Y,ETA,XX,AMU)*( AMU-A)/(AMU* 
1TANH( AMU*H)-ANU) 

Q2X(Y , ETA,X,XI,AMU)=-P2X(Y,ETA,X,XI,AMU)*( AMU-A) /(AMU* 
1TANH( AMU*H)-ANJ) 

Q2Y(Y,ETA,X,XI,AMU)=-P2Y(Y,ETA,X,XI,AMU)*( AMU-A) /(AMU* 
1TANH( AMU*H)-ANU) 

Q2YY( Y,ETA,X,XI, AMU)=-P2YY (Y,ETA,X,XI , AMU ) * ( AMU- A ) / 
1( AMU*TANH( AMU*H) -ANU) 

Q10(Y ,ETA,XX)=-COSH(A*(Y+H) ) *COSH ( A* ( ET A+H) )*COS(A*XX) 
1*A/(C0SH( A*H)**2*( ( A*A-ANU*ANU)*H+ANU)) 

Q2CX(Y,ETA,X,XI ) =A*COSH( A*( ETA  +  H ) ) *S IN( A* ( X-X I ) )*COSH 
1 ( A*( Y+H ) )*A/i:OSH( A*H)**2*( (A* A- ANU* AN J )*H+ANU) ) 

Q20Y( Y,ETA,X,XI ) =-A*COSH( A* ( ETA+H ) ) *S INH ( A* ( Y+H) ) *COS 
l(A*(X-XI))*A/(C0SH(A*H)**2*(l A* A- ANU* ANU ) *H+ANU) ) 

Q20YY(Y,ETA,X,XI ) =-A*A*COSH ( A* ( Y+H) )*COSH(A*( ETA+H))* 
1C0S(A*( X-XI ) )*A/(COSH( A*H)**2*( ( A* A- ANU* ANU ) *H+ANU ) ) 

Q1S(Y,ETA,XX,AMU)=-P1 ( Y , ETA , XX , AMU) / ( A0+ AMU*H*( AA +BB* 
1 (AMU-A) +CC*( AMU-A) ** 2 +DD*( AMU-A )**3) ) 

Q2XS(Y, ETA,X,XI, AMU ) =-P2X ( Y , ETA , X ,XI , AMU )/( AO+AMU*H* 
1 (AA+BB* (AMU-A) +CC*( AMU-A ) **2+QD* ( AMU-A ) **3) ) 

Q2YS( Y,ETA,X,XI , AMU) =  -P2Y ( Y ,ETA , X ,XI , AMU) / ( AO+AMU*H* 
1 (AA+BB* ( AMU-A )+CC*( AMU-A ) **2+DD* ( AMU-A ) **3 ) ) 

Q2YYS(Y,ETA,X,XI  , AMU ) =-P2YY ( Y , ETA, X, X  I , AMU ) / ( A0+ AMU*H* 
1( AA+BB* ( AMU-A )+:C*( AMU-A )**2+DD*( AMU-A )**3) ) 

FUNKY,  ETA,XX,AMU)=EXP(-AMU*H)*SINH(AMU*ETA)*SINH 
1( AMU*Y)*COS( AMU*XX)/( AMU* COSH ( AMU*H) ) 

FUN2( Y, ETA,XX)=4.*3.14159*C0SH( A* (Y+H) )*CCSH( A* 
1( ETA+H) )*COS( A*XX)/( 2. *A*H+ SH2AH) 

FUNR(X, Y,XI,ETA,ANX,ANY)=( (X-XI ) *ANX+( Y+ET A ) *ANY ) / 
1( (X-XI)**2+(Y+ETA)**2) 

FUN3X(Y,ETA,X,XI , AMU ) =EXP ( -AMU*H ) *S INH( AMU* ET A) *S I NH 
11 AMU*Y)*SIN( AMU* (X-X I) ) /COSH ( AMU *H ) 

FUN3Y(Y,ETA,X,XI , AMU) =-EXP ( -AMU*H) *SINH ( AMU*ETA ) *COSH 
1 ( AMU*Y)*COS( AMU* (X-XI ) ) /COSH(  AMU*H ) 

FUN3YY( Y,ETA,X,XI,AMU)=-EXP(-AMU*H)*SINH( AMU*ETA)* 
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1SINH( AMU*Y)*AMU*AMU*COS(AMU*<X-XI ) ) /COSH ( AMU*H ) 

FUN4(Y,ETA,XX,ANX,ANY,SIGN)=4.*3.  14159*  A*COSH  (  A* 
i(  ETA+H)  )*(  ANY*SINH(  A*(Y+H)  )  *CD S  (  A*XX)  -ANX*C OSH(  A* 
2(Y+H))*SIN(A*XX)*SIGN)/(2.*A*H+SH2AH) 

EVALUATION    OF    THE    FINITE    INTEGRAL    IN    G    TO    DETERMINE 
THE    SIZE    OF    THE    SUBDIVISIONS 


10 


30 

40 


45 


50 

60 


XX=ABS(X-X 
IF(X.LT.XI 
IFiX.GT.XI 
IFU.EQ.XI 
DO  50  N  =  l, 
DELMU=2*A/ 
AMU=0.0 
SUM=0.0 
F0=(Q1( Y,E 
LL=6*N+3 
DO  40  NN=1 
IF( ABS( AMU 
IFIA3S( AMU 
IFUBSiAMU 
IF(ABS( AMU 
F1=(Q1(Y,E 

1DELMU/3.-A 
F2=(Q1 (Y,E 

l(AMU+2.*DE 
F3=(Q1( Y,E 

1-A) 
GO    TO    30 
F0=(Q1( Y 

1-A) 
GO    TO    40 
SUM=  (DE 

F0=F3 

AMU=AMU+DE 
TEST(N)=SU 
IF(N-l)     50 
MN=N-1 
MM=6*MN+3 
IF(ABS( (TE 
CONTINUE 
CONTINUE 
PV1F=2.*SU 


I) 

)    SI3N=-1.0 

)    SIGN=1.0 

)     SIGN=0.0 

15 

(6*N+3) 


TA,XX, AMU)-Q10(Y,ETA,XX) )/< AMU-A) 

,LL 

-A)  .LT. .00001)     GO    TO     10 

+DELMU-A) .LT.. 00001)     GO    TO    10 

+DELMU/3.-A).LT. .00001  )    GO    TO     10 

+2.*DELMU/3.-A) .LT..00001)     G3    TO    10 

TA,XX, AMU+DELMU/ 3. )-Q10(Y,ETA,XX)  )/(AMU+ 

) 

TA,XX, AMU+2.*DELMU/3.  ) -Q10 ( Y , ETA, XX ) ) / 

LMU/3.-A) 

TA,XX, AMU  +  DELMU )-Q10(Y,  ETA, XX) )/( AMU+DELMU 


ETA, XX, AMU+DELMU J -Q10( Y, ETA, XX) ) /( AMU+DELMU 


LMU/3. )*(F0+3.*F1+3.*F2+F3)+SUM 
LMU 

M 

',50,45 


ST(N)-TEST(N-1) ) /TEST(N) ) .LT. .010)  GO  TO  60 


EVALUATION  OF  THE  FINITE 
2*A/MM  SUBDIVISION  SIZE 


INTEGRAL  IN  GN  USING 


DELMU=2*A/MM 

AMU=0.0 

SUMX=0.0 

SUMY=0.0 

F0=(Q2X(Y,ETA,X,XI ,AMU)-Q2  0X(Y ,ETA,X,XI ) )/{ AMU-A) 

Y0=(Q2Y(Y,ETA,X,XI , AMU) -Q20Y ( Y , ETA ,X , XI ) ) /( AMU-A) 

DO  80  NN=1,MM 

IF(ABS( AMU-A) .LT. .00001)  GO  TO  70 

IF(ABS( AMU+DELMU-A) .LT.. 00001)  GO  TO  70 

IF(ABS( AMU+DELMU/3.-A) .LT..00001 )  GO  TO  70 

IF(ABS( AMU  +  Z.OELMU/3.-A) .LT. .00001)  GO  TO  70 

F1  =  (Q2X(Y, ET A, X, XI, AMU+DELMU/ 3. ) -Q20X ( Y , ETA, X , X  I ) )/ 
1{ AMU+DELMU/3.-A) 

F2=(Q2X(Y,ETA,X,XI , AMU+2.*DE LMU/3 .) -Q20X ( Y , ETA , X , XI ) )/ 
1( AMU+2.*DELMU/3.-A) 

F3=(Q2X(Y,ETA,X,XI, AMU  +  DE LMU  )  - Q2  OX  { Y  ,  ET  A  ,  X  ,  X  I  )  )  / 
1< AMU+DELMU-A) 

Y1=(Q2Y(Y,ETA,X,XI , AMU+DELMU/ 3 . )-Q20Y(Y ,ETA,X,XI) )/ 
1( AMU+DELMU/3.-A) 

Y2=(Q2Y(Y,ETA,X,XI,AMU+2.*DELMU/3.)-Q20Y(Y,ETA,X,XI) )/ 
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l(AMU+2.*DELMU/3.-A) 

Y3={Q2Y(Y,ETA,X,XI,AMU+DELMU)^Q2  0Y(Y»ETA,X,XI ))/ 
1( AMU+DELMU-A) 

GO  TO  75 
70  CONTINUE 

F0=(Q2X(Y ,ETA,XtXI , AMU+DELMU ) -Q2CX( Y , ET A , X , XI ) )/ 
1( AMU+DELMU-A) 

Y0=(Q2Y(Y,ETA,X,XI,AMU+DELMU)-Q2GY(Y,ETA»X,XI  )  )/ 
1 (AMU+DELMU-A) 

GO    TO    80 
75    CONTINUE 

SUMX=SUMX+  (DELMU/8.  )  *(  FO  +  3  .*F  1  +  3  .*F2+F3  J 

SUMY=SUMY+  (DELMU/8. )* ( YO+3 . *Y 1+3.*Y2 +Y3 ) 

F0  =  F3 

Y0  =  Y3 
30  AMU=AMU+DELMU 

PV2FX=2.*SUMX 

PV2FY=2.*SUMY 

EVALUATION  OF  THE  FINITE  INTEGRAL  IN  GYY  USING 
2*A/MM  SUBDIVISION  SIZE 

DELMU=2*A/MM 

AMU=0.0 

SUMYY=0.0 

YY0=(Q2YY(Y,ETA,X,XI  ,AMU) -Q20Y Y ( Y  ,ETA , X , X  I )  )/( AMU-A) 

DO    350    NN=1,MM 

IF(ABS( AMU-AJ .LT. 0.00001)     GO    TO    320 

IF     (A3S( AMU+DELMU-A) .LT. 0.00001 )    GO    TO    320 

IF    UBSCAMU+DELMU/3. -A). LT. 0.00001)    GO    TO    320 

IF    (ABS(AMU+2.*0ELMU/3. -A) .LT. 0.00001)     GO    TO    320 

YY 1= ( Q2 YY ( Y, ET A, X, XI T AMU+DELMU /3.)-Q20YY(Y, ETA, X, XI ) )/ 
1 ( AMU+DELMU/3. -A) 

YY2=(Q2YY( Y,ETA,X,XI , AMU+2 . *DE LMU/3. ) -3  20YY (Y,ETA,X, 
IX  I) )/( AMU+2.*DELMU/3.-A) 

YY3=(Q2YY(Y,ETA,X,XI , AMU+DELMU ) -Q2QYY ( / , ETA,X,XI ) )/ 
1 (AMU+DELMU-A) 

GO  TO  325 
320  CONTINUE 

YYO=(Q2YY(Y,ETA,X,XI,AMU+DELMU)-Q20YY(Y,ETA,X,XI) )/ 
1( AMU+DELMU-A) 

GO  TO  350 
325  CONTINUE 

SUMYY=SUMYY+( DELMU/8. )*( Y YO+3. *YY 1+3 .*YY 2+YY3 ) 

YY0=YY3 
350  AMU=AMU+DELMU 

PV2FYY=2.*SUMYY 

EVALUATION  OF  THE  INFINITE  INTEGRAL  IN  G,  GNV  AND  GYY 
SIMULTANEOUSLY 

AMU=2*A 

DELMU0=0ELMU 

F0=Ql(Y,ETA,XX,AMU)/( AMU- A) 

F0X=Q2X (Y,ETA,X,XI , AMU) / ( AMU-A ) 

F0Y=Q2Y(Y,ETA,X,XI ,AMU)/( AMU-A) 

F0YY=Q2YY( Y,ETA,X,XI ,AMU) /( AMU-A) 

SUM=0.0 

SUMX=0.0 

SUMY=0.0 

SUMYY=0.0 

DO  100  NN=1,200 

DO  95  N  =  l, 20 

F1=Q1(Y, ETA, XX, AMU+DELMU/3. ) / ( AMU+DELMJ/ 3 .- A) 

F1X=Q2X(Y,ETA,X,XI , AMU+DELMU/3 . ) /( AMU+DELMU/3 .-A) 

F1Y=Q2Y(Y,ETA,X,XI , AMU+DELMU/3 . )/( AMU+DELMU/3. -A) 

F1YY=Q2YY( Y, ETA, X, XI, AMU+DELMU/3. ) / ( AMJ+DELMU/3 .-A ) 

F2=Ql(Y,ETA,XX,AMU+2.*DELMU/3. ) / ( AMU+2. -DELMU/3. -A ) 
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F2X=Q2X(Y,ETA,XTXI,AMU+2.*DELMU/3. ) / ( AMU* 2. *DELMU/3 . 

1-A) 
F2Y=Q2Y(Y,ETA,X,XI , AMU*2.*DE LMU/3 . ) / ( AMU* 2 . *DELMU/3 . 

1-A) 
F2YY=Q2YY{ Y,ETA,X,XI,AMU*2.*DELMU/3.)/( AMU*2.*DELMU/3. 

1-A) 
F3=Q1(Y,ETA,XX,AMU*         DELMU ) / { AMU+DELMJ- A ) 
F3X=Q2X(Y,ETA,X,XI , AMU+DELMU) / ( AMU+DELMU-A) 
F3Y=Q2Y<Y,ETA,X,XI,AMU+DELMU)/( AMU  +  DELVIU- A  ) 
F3YY=Q2YY( Y,ETA,X,XI, AMU*DELMU)/( AMU*DELMU-A) 
SUM=  (OELMU/3. )*(F0*3.*F1+3.*F2*F3)+SUM 

SUMX=  (DELMU/8.)*(F0X*3.*F1X*3.*F2X*F3X)*SUMX 

SUMY=  <DELMU/8.)=MF0Y  +  3.*F1Y  +  3.*F2Y*F3Y)*SUMY 

SUMYY=SUMYY*(0ELMU/8. ) * ( FOYY+3 .* F1YY*3 . *F2Y Y+F3YY ) 
FO  =  F3 
F0X=F3X 
FOY=F3Y 
F0YY=F3YY 

ASM=EXPUMU*(ETA*Y) )/AMU 
IF(ASM.LT..OOOi )     GO    TO    86 

95  AMU-AMU+OELMU 
86    T£ST(NN)=SUM 

TESTT(NN)-SUMX 
TESTTT(NN)=SUMY 
TESTYY(NN)=SUMYY 
IF(NN-l)     100,100,97 

97  CONTINUE 

IFUSM.LT.  .00010)     GO    TO    105 
IF(A6S(SUM) .LT. .0001)     GO    TO    98 

IF (ABS ( (TEST(NN) -TEST (NN-1 )  ) /TEST  INN)) .GT..001) 
1G0    TO    100 

98  IF(ABS( SUMX) .LT..0001)     GO    TO    99 

IF  UBS  < (TESTT(NN)-TESTT(NN-n ) / TESTT(NN) ) .GT. .001) 
1G0    TO    100 

99  IF(ABS( SUMY) .LT. .0001)     GO    TO    93 

IF(ABS( (TESTTT(NN)-TESTTT(NN-l) ) / TESTTT ( NN) ) . GT. . 001 ) 
1G0    TO    100 
93    IF(ABS( SUMYYJ.LT. 0.0001)     GO    TO    96 

IF    (ABS( { TESTY Y (NN)-TESTYY*  NN-1)  ) /TESTY Y ( NN  ) )  .GT..001) 
1G0    TO    100 

96  CONTINUE 
GO    TO    105 

100    CONTINUE 

102  WRITE(6,103) 

103  F0RMAK3X34HINFINITE    INTEGRAL    DID    NOT    CONVERGE) 
105    CONTINUE 

PV1I=2.*SUM 

PV2IX  =  2  .*SUMX 

PV2IY=2.*SUMY 

PV2IYY=2.*SUMYY 

DELMU=.2/H 

IF  (J.EQ.25)  GO  TO  240 

AMU=0.0 

SUM=0.0 

SUMX=0.0 

SUMY=0.0 

SUMYY=0.0 

F0=0.0 

F0X=FUN3X{ Y,ETA,X,XI,AMU) 

F0Y=FUN3Y{ Y,ETA,X,XI ,AMU) 

F0YY=FUN3YY( Y, ETA,X,XI ,AMU) 

DO  200  NN=1,100 

IF( AMU*H.GT.5. )  DELMU=.l 

DO  195  N=l,20 

F1X  =  FUN3X( Y,ETA,X,XI , AMU*DE LMU/3.  J 

F2X=FUN3X( Y,ETA,X,XI , AMU* 2. *DE LMU/3. ) 

F3X  =  FUN3X( Y,ETA,X,XI , AMU  +  DELMU  ) 

F1Y  =  FUN3Y{ Y,ETA,X,XI , AMU  +  DE LMU/3  .  ) 

F2Y=FUM3Y( Y , ETA , X , XI , AMU+2. *DE LMU/3. ) 

F3Y=FUN3Y( Y,ETA,X,XI , AMU+DELMU) 

F1YY  =  FUN3YY(Y,ETA,X,XI , AMU+DELMU/ 3  . ) 

F2YY=FUN3YY(Y,ETA,X,XI , AMU+2 .*DE LMU/3 . ) 
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F3YY=FUN3YY< Y , ETA ,X ,X I , AMU+DELMU ) 

IF< AMU+DELMU. LT. .001)     GO    TO    120 

F1=FUN1 (Y, ETA, XX, AMU+DELMU/ 3.) 

F2=FUNi(Y,ETA,XX,AMU+2.*DELMU/3.) 

F3=FUN1(Y, ETA, XX, AMU+DELMU) 

GO  TO  130 
120  CONTINUE 

Fl=EXP(-<  AMU  +  DELMU/ 3.  )*H)*COS(  ( AMU  +  DELMU /3.  )*XX)* 
1( AMU+DELMU/3. )*Y*E T A/COSH  I ( AMU  +  DELMU/ 3. )*H) 

F2=EXP(-{ AMU*2.*DELMU/3.)*H)*CGS(  ( AMU+2 . *DELMU/3 .  )*XX) 
l=MAMU  +  2.*DELMU/3. ) *Y* ETA/ COSH ( ( AMU+2 .*DE LMU/3 . ) *H ) 

F3=EXP(-( AMU+DELMU )*H)*COS( ( AMU+DELMU )* XX )*( AMU+DELMU) 
1*Y*ETA/C0SH( (AMU+DELMU)*H) 
130  CONTINUE 

SUM=    (OELMU/8. ) * ( FO +3 .*F1 +3 .*F2 +F3 ) +SUM 

SUMX=  (DELMU/8.)*( FOX+3 . *FlX  +  3  .*F2X  +  F3X ) +  SUMX 

SUMY=  (DELMU/8.)*(F0Y+3.*F1Y+3.*F2Y+F3Y)+SUMY 

SUMYY=SUMYY+(DELMU/8. ) *( FOY Y+3 . *F 1YY+3 . *F2YY+F3YY ) 

F0=F3 

F0X=F3X 

F0Y=F3Y 

F0YY=F3YY 

ASM=EXP(AMU*(ETA+Y) ) 

IF(ASM.LT..0001)  GO  TO  196 

195  AMU=AMU+DELMU 

196  TEST(NN)=SUM 
TESTT(.NN)=SUMX 
TESTTT( NN)=SUMY 
TESTYY(NN)=SUMYY 
IF(NN-l)  200,200,199 

199  CONTINUE 

IFtASM.LT.  .00010)  GO  TO  205 
IF(ABS( SUM) .LT.. 0001)  GO  TO  2C6 

IF(ABS( ( TEST(NN)-TEST(NN-1) ) /TEST ( NN) ) .GT. . 0010) 
1G0  TO  200 

206  IF(ABS( SUMX) .LT. .0001)  GO  TO  207 

IFUBSt  ( TESTT(NN)-TESTT(NN-1) )/TESTT(NM) ) .GT..0010) 
1G0  TO  200 

207  IF (ABS( SUMY) .LT. .0001)  GO  TO  204 

IF (ABS <  ( TESTTT(MN)-TESTTTdMN-l)  )/TESTTT(NN)  ).GT. .0010) 
1G0  TO  200 

204  IF(ABS( SUMYY) .LT. 0.0001)     GO    TO    208 

IF    (ABS(  (  TESTY  Y  (  iMN  )-TE  STY  Y(NN-l)  )  /  UESTY  Y  (  NN)  )  )  . 
1GT..001 )    GO    TO    200 

208  CONTINUE 
GO    TO    205 

200  CONTINUE 
WRITE{6,202) 

202    FORMAT (3X14HN0    CONVERGENCE) 

205  CONTINUE 
GINF=-2*SUM 
GXINF=2.*SUMX 
GYINF=2.*SUMY 
GYYINF=2.*SUMYY 
GNSING= .5 

IF    (I.EQ.25)    GO    TO    218 

IF(I.EQ.J)    GO    TO    220 

AIJJ-I-J 

THETA1=ABS( AIJJ)*DELTHE 

AINJJ=I+NPTS-J 

THETA2=ABS(AINJJ )*DELTHE 

AJNII=I-J-NPTS 

THETA3=ABS( AJNI I )*DELTHE 

THETA=THETA1 

IF1THETA2.LT.THETA)  THETA=THETA2 

IF(THETA3.LT.THETA)  THETA=THETA3 

IF(THETA.GT..15)  GO  TO  218 

GS I NG=D ELT HE *ALOG( 2. ) +  2  .  *( -DEL  THE/2.  +  OELTHE/ 4.  + 
1THETA/2. )*ALQGi  THETA/ 2 . +DELTHE/4 . )-( THET A/ 2 .-DEL THE/ 4. 
2) *ALOG( THET A/2. -DELTHE/4. ) - ( DELTHE/4 .+THET A/ 2. i ** 3/ 18. 
3+ (THET  A/ 2. -DEL  THEM. )**3/18.-( THETA/2. +DELTHE/+ . )**5/ 
4( 180.*5.)+(THETA/2.-DELTHE/4. )**5/< 180.*5. ) -(THET A/2. + 
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5DELTHE/4.)**7/( 2  83  5.*7.)+( THE T A/ 2. -DELTHE/4 . ) **7/ 
6(2835. *7.)  ) 

GXSING=(X-XI )/( ( X-XI )**2+( Y-ETA )**2) 

GYSING=( Y-ETA) /( (X-XI )**2+( Y-ETA )**2) 

GO  TO  230 
218  GSING=DELTHE*ALOG(SQRT( ( X-X I ) **2 ♦ ( Y-ETA ) **2 ) ) 

GXSING=(X-XI )/( (X-XI )**2+( Y-ETA )**2) 

GYSING=(Y-ETA)/( (X-XI ) **2+( Y-ETA ) **2 ) 

GO    TO    230 
220    CONTINUE 

GSING=DELTHE*AL0G(DELTHE/2. ) -DELTHE- ( DE LTHE**3) / 
1(  l8.*16.)-(DELTHE**5)/( 180 . * 5. *256. ) -  (  3ELTH E**7 ) / 
2(2835. *7.*256.*16.  ) 

GXSING=GNSING*ANX 

GYSING=GNSING*AMY 
230    CONTINUE 

G      =GSI NG/DELTHE-ALOG(SQRT( ( X-X I ) **2  + ( Y +ET A  )**2 ) )+PVlF 
H-PV1H-GINF+CMPLX(0.,-1.  )*FUN2( Y  ,  ETA, XX) 

GX=GXSING-FUNR(X, Y,XI  ,ETA,1  .0  , 0  .0  )  +PV2FX  *PV2IX*-GXI  NF  + 
1FUN4(Y ,ETA,XX,1 . 0 ,0.0, SIGN )*CMPLX (0.0,-1 .0) 

GY=GYSING-FUNR( X,Y,XI ,  ETA,0 .0,1 .0 ) +P V2F Y +  PV 21 Y  +  GY I NF+ 
IF UN4(Y, ETA, XX, 0.0, 1.0, SIGN) *CMPLX (0.0,-1.0) 

IF    (I.LE.NPTS)     30    TO    235 

FURYY=1./SQRT( ( (X-XI ) **2*( Y-ET A ) **2 ) **3 ) - ( ( Y-ETA)** 2+ 
1( Y-ETA) )/SQRT( ( (X-XI » **2+( Y-ETA) **2)**5 ) 

FUNYY=1./SQRT(  (  (X-XI  )  **2  +  (  Y  +  ETA  )  **2)  **3  )  -  (  (  Y+ETA  )  **2  + 
KY  +  ETA)  )/SQRT(  (  (X-XI  )  **2*(  Y  +ET  A  )  **2)  **5  ) 

G YY=FURYY+FUNY Y +PV2FY Y+PV2 I YY+GYY INF *CMPLX( 0.0,-1.0)* 
1FUN2(Y ,ETA,XX)*A*A 
235  CONTINUE 

IF  (I.EQ.25)  GO  TO  250 

GN=GNSING-FJNR(X,Y,XI , ETA , A NX , ANY ) +{ PV2F X  +  P V2 IX+GX INF ) 
l*ANX+( PV2FY+PV2IY+GYINF)*ANY+FUN4(Y,ETA,XX, ANX,ANY, 
2SIGN)*CMPLX(0. 0,-1.0) 

GO  TO  250 
240  CONTINUE 

G=-(PV1F4-PV1I)+CHPLX(0.  ,-l.)*FUN2(Y,ETA,XX) 

GN=-( PV2FX+PV2IX)*ANX-(PV2FY+PV2IY)*ANY*FUN4( Y, ETA, XX, 
1 ANX, ANY , SI GN)*CMPLX(0. 0,-1.0) 
250  CONTINUE 

RETURN 

END 

THIS  SUBROUTINE  CALCULATES  G  AND  ITS  DERIVATIVES, 
GX,GY,GN,  AND  3YY,  BY  USE  OF    THE  SERIES  FORM  FOR  THE 
CASE  (X(I)  -  X(J))  GREATER  THAN  SMIN 

SUBROUTINE  GREENS  ( A , ANU, SH2 AH, SHAH ,CHAH ,COSAMU , S INAMU 
1, AMU,COEFG,SHYI , CHY I , SHY  J , C HY J , I  , J , XV , Y V , XC , YC , AN  XI , 
2ANYI,  ANXJ,A,MYJ,GIJ,GNIJ,G.MJI  ,  oX  Ij  ,GX  J  I  ,  G  YI  J  ,GY  J I  ,GYY) 

COMPLEX  GIJ,GNIJ,GNJI ,GXIJ,GXJ I ,GYIJ,GYJI , GYY 

DIMENSION  COS AMU ( 200,2  5) ,SINAMU( 200,25)  , AMU (200)  , 
1C0EFG(200) 

DIMENSION  TESTK40)  ,TEST2(40)  ,TEST3(40)  ,TEST30(40)  , 
1TEST4(40) 

COMMON/ CON ST /H,D,DELTHE,S MI N,NPTS,NSPTS 

CCMMON/CPI/PI 

SUM1=0.0 

SUM2=0.0 

SUM3=0. 0 

SUM30=0.0 

SUM4=0.0 

DO  20  N=l,40 

DO  15  KK=1,5 

K=(N-1) *5  +  KK 

SNI=SINAMU(K,I ) 

SNJ-SINAMU(K,J) 

CSI=COSAMU(K,I) 

CSJ=COSAMU(K,J) 
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AK=AMU(K) 

CCK=COEFG(KJ 

VAL=AK*ABS<XV-XC) 

IF    (VAL.GE.75.0)     GO    TO    10 

EXIJ  =  EXP(-VAU 

GO    TO    11 

10  EXIJ=0.0 

11  CONTINUE 
SUM1=SUMH-C0K*CSI*CSJ*EXIJ 
SUM2=SUM2+CCK*CSI*CSJ*EXIJ*AK 
SUM3=SUM3+C0K*AK*SNI*CSJ*EXIJ 
SUM3OSUM30+C0K*AK*SNJ*CSI*EXIJ 
SUM4=SUM4+AK*AK*C0K*CSI*CSJ*EXIJ 

15  CONTINUE 
TEST1(N)=SUM1 
TEST2(N)=SUM2 
TEST3(N)-SUM3 
TEST30(N)=SUM30 
TEST4(N)=SUM4 
IF(N-l)     20,20,16 

16  IF(ABS( SUM1) .LE. 0.000001)     30    TO    17 

IF  UBS  I  (TESTl(.M)-TESTUN-l)  J/TESTKN)  ).GT..0010)  GO 
1T0  20 

17  IF(ABS< SUM30) .LE .0.000001)  GO  TO  18 

IF<ABS(  (TEST30(N)-TEST30(N-1)  )/TEST30(N)  )  .GT.  .0010) 
1G0  TO  2  0 

18  IF(ABS( SUM3) .LE. 0.000001)  GO  TO  19 
IF(ABS((TEST3(N)-TEST3(N-1))/TEST3(N)).GT..0010)  GO 

1T0  20 

19  IF(A6S< SUM2) .LE. 0.000001)  GO  TO  21 

IF(ABS( ( TEST2(N)-TEST2(N-1) )/TEST2(N)) .GT..0010)  GO 
1T0  20 
21  IF(ABS( SUM4) .LE. 0.000001)  GO  TO  30 

IF(  ABS(  (  TEST4<N)-TEST4(N-1  )  )/TEST4(N)  )  .GT..0010)  GO 
1T0  20 

GO  TO  30 

20  CONTINUE 
WRITE(6,25)  I, J 

25  FORMAT(10X23HGREENS  DID  NOT  CONVERGE , 2X2HI  =  I  2 , 

12X2HJ=I2) 
30  CONTINUE 

IF  (XV.LT.XC)  SIGN=-1.0 

IF  (XV.GT.XC)  SIGN=1.0 

IF  (XV.EQ.XC)  SIGN=0.0 

TERM4=4.*PI*CHYI*CHYJ*SIN( A*ABS(XV-XC)) / < 2. *A*H+SH2AH) 

TERM5=4.*PI*CHYI*CHYJ*CCS( A*ABS( XV-XC) )/ ( 2 . *A*H+SH2AH) 

TERM6=A*TERM5*SIGN 

TERM7=A*TERM4*SIGN 

TERM8=4.*PI*A*SHYI*CHYJ*C0S(A*A3S(XV-XC )  )/<2.*A*H  + 
1SH2AH) 

TERM9=4.*PI*A*SHYI*CHYJ*SIN(A*ABS(XV-XC) )/( 2.*A*H+ 
1SH2AH) 

TERM80  =  4.*PI*A*SHYJ*CHYI*C0S(A*ABS(X\/-XC ) )/(2.*A*H+ 
1SH2AH) 

TERM90=4.*PI*A*SHYJ*CHYI*SIN(A*ABS(XV-XC) )/(2.*A*H+ 
1SH2AH) 

TERM10=A*A*TERM4 

TERM11=A*A*TERM5 

IF    (J.EQ.25)     GO    TO    70 

GIJ=CMPLX( SUM1+TERM4,-TERM5) 

GXIJ=    CMPLX(-SIGN*SUM2*TERM6,TERM7) 

GYIJ=CMPLXl-SUM3+TERM9,-TERM8) 

IF    ( I.EQ.25)     GO    TO    40 

GXJI=CMPLX(SIGN*SUM2-TERM6,-TERM7) 

GYJI=CMPLX(-SUM30*TERM90,-TERM80) 

GNI J=CMPLX(-ANXI*SIGN*SUM2+ANXI*TERrt6-ANYI*SUM3*ANYI* 
lTERM9f ANXI*TERM7-ANYI*TERM8) 

GNJI=CMPLX(ANXJ*SIGN*SUM2-ANXJ*TERM6-ANYJ*SUM3G+ANYJ* 
1TERM90,-ANXJ*TERM7-ANYJ*TERM80) 

GO    TO    60 
40    CONTINUE 

GYY=CMPLX(-SJM<H-TERM10,-TERM11 ) 
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60  CONTINUE 

GO  TO  80 
70  GIJ=CMPLX(-SUM1-TERM4,-TERM5) 

GNIJ=CMPLX(ANXI*SIGN*SUM2-ANXI*TERM6+ANYI*SUM3-ANYI' 
LTERM9fANXI*TERM7-ANYI*TERM8) 
80  CONTINUE 

RETURN 

END 

THIS  SUBROUTINE  INVERTS  COMPLEX  MATRICES  TO  SOLVE 
THE  MATRIX  EQUATION  ALPHA( I . J } *F { I , 1)  =  HH(I,1)  AND 
ALPHAU  ,J)*F1(  I  ,1)  =  PK(I,1) 

SUBROUTINE  COMAT ( N  ,M , A , B , 0 , I ) 

INTEGER  C»H,R,Q,Z 

COMPLEX  A,B,D,TT,P 

DIMENSION  A(NtN) ,6(N,M) ,C( 100,3) 

D  =  ( 1 .0,0.0) 

DO  20  J  =  1,N 

20  C(J,3)  =  0 

DO  21  K  =  1,N 

TT  =  (0.0,0.0) 

T  =  0.0 

DO  4  J  =  1,N 

IF  (C( J, 3)  .EQ.  1)  GO  TO  4 

DO  5  H  =  1,N 

IF  (C(H,3)  -  1)  15,5,12 

15  IF  (T  .GE.       CABS( A( J,H) ))   GO  TO  5 
R  =  J 

Q  =  H 

T  =       CABS( A( J,H) ) 
5  CONTINUE 
4  CONTINUE 

C(Q,3)  =  C(Q,3)  ♦  I 

C(K,1)  =  R 

C(K,2)  =  Q 

IF  (R  . EQ.  Q)  GO  TO  11 

D  =  -D 

DO  8  L  =  1  ,N 

TT  =  AIR  ,L) 

A(R,L)  =  A(Q,L) 
8  A(Q,L)  =  TT 

IF(M  .LE.  0)  GO  TO  11 

DO  2  L  =  1  ,M 

TT  =  B( R,L) 

B(R,L)  =  B(Q,L) 

2  BIQ,L)  =  TT 
11  P  =  A(Q,Q) 

A(Q,Q)  =  (1.0,0.0) 
DO  13  L  =  1,N 
13  A(Q,L)  =  A(Q,L)/P 

IF  (M  .LE.  0)  30  TO  1 
DO  3  L  =  1,M 

3  B(Q,L)  =  B(Q,L)/P 
1  DO  21  I    =  1,N 

IF  (Z  .EQ.  Q)  GO  TO  21 
TT  =  A( Z,Q) 
A(Z,Q)=(0. 0,0.0) 
DO  16  L  =  1,N 

16  A(Z,L)  =  A(Z,L)  -  A<Q,L)*TT 
IF  (M  . LE.  0)  GO  TO  21 

DO  17  L  =  1,M 

17  B(Z,L)  =  B(Z,L)  -  B(Q,L)*TT 

21  CONTINUE 

DO  19  II  =  1,N 

L  *  N  +  1  -  II 

IF(C(L,1)  .EQ.  C(L,2))  GO  TO  19 

R  =  C(L,1) 

Q  =  C(L,2) 
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DO  7  K  *  l. N 

TT  =  A(K,R) 

A(K,R)  =  A(K,Q) 
7  A(KfQ)  =  TT 
19  CONTINUE 

DO  18  K  =  UN 

IF  (C(K,3)  .NE.  1)  GO  TO  12 
18  CONTINUE 

I  =  1 
50  RETURN 
12  I  =  2 

GO  TO  50 

END 


//GO.SYSIN    DD  * 
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